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INTRODUCTION AND SUMMARY 


7 * 

This report describes the SINFLO modification package for SINDA, which 
was developed by the Vought Systems Division (VSD) of LTV Aerospace Corporation 
under subcontract to TRW Systems Group during the period of August 1973 to 
February 1975. Also included in this report is a description of the capabilities 
added during the development of SINDA-VERSION 9^ . The SINFLO package was developed 
to modify the SINDA preprocessor to accept and store the input data for fluid 
flow systems analysis and adding the FLOSOL user subroutine to perform the flow 
solution. This greatly reduced and simplified the user input required for analysis 
of flow problems. Also, a temperature calculation method, the Flow-Hybrid method 
which was developed in previous VSD thermal simulator routines , was incorporated 
for calculating fluid temperatures. The calculation method accuracy was improved 
by using fluid enthalpy rather than specific heat for the convective term of the 
fluid temperature equation. 

The effort described herein was performed under Task Order LT-1 and Task 
Order LT-2 of Subcontract 183LK3E of NASA Contract NAS9-10435. Task Order LT-1 
calls for the completion of the following tasks: 

A. Optimize the SINDA Routine flow input data which includes the 
following effort: 

1. Establish the best input format for the flow systems data 
block to be added. 

2. Establish the format for storing the flow systems data by 
the preprocessor. 

3. Submit the user input and preprocessor output formats for 
mutual agreement between NASA-JSC, TRW, Inc. and VSD. 

4. Modify the preprocessor to accept the input from the new 
flow data block and store for use with processor routines. 

B. Write routine to perform fluid flow analysis using data stored 
by proprocessor. 

C. Develop Fluid Hybrid Routine for SINDA which will take advantage 
of the temperature calculation equation form for the fluid lumps 
to calculate the fluid and the tube lumps temperatures using an 
implicit method while the remaining structure lumps will be 
calculated utilizing the method specified by the user. 


*Superscripts indicate reference numbers in Section 7.0 
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Task Order LT-2 required the completion of the following tasks: 

A. Modify the following SINDA execution routines to interface 
with subroutine FLUID so that thermal analysis of fluid 
flow systems may be performed: (1) SNFRDL, (2) SNFRWD, 

(3) CINDSL, (4) FWDBCK, (5) SND5NR, (6) STDSTL. 

B. Add the capability to analyze valves which will split 
incoming flow between the two outlet sides of the valve 
in proportion to the valve position regardless of the 
pressure balance. 

These tasks were completed and the resulting routines added to and 
substituted into the SINDA general thermal analyzer routine expanded the capa- 
bilities of the SINDA to include analysis of systems containing flowing fluids, 
fluid system controls and heat exchangers. A pressure-flow analysis of a system 
containing an arbitrary tube network is performed simultaneously with the thermal 
analysis during transient or steady state solutions. This permits the mutual 
influences of thermal and fluid problems to be included in the analysis. 

The general flow solution capabilities include extensive valve character- 
izations and ability to match pump curves and system presrure-flow characteristics. 
The valves have been formulated so that either cooling (space radiator) or heating 
(solar absorber) situations may be controlled with any of the valve types. Pump 
options included are pressure rise as a tabulated function of system flow rate and 
pressure rise as a polynomial function of flow rate. 

The formulation of the capabilities added during this effort are described 
in Section 2.0, modifications to the SINDA subroutines are described in Section 3.0, 
and the data input requirements for the new data block are described in Section 4.0. 
Section 5.0 describes user subroutines which have been added or modified by VSD 
including those developed for SINDA-VERSI0N9. 

Appendix A contains a description of the capabilities incorporated into 
subroutines during the development of SINDA-VERSION 9 to facilitate analysis of 
radiation heat transfer in an enclosure. A description of the usage of the plot 
program is presented in Appendix C. 

A discussion of the flow data storage is presented in Appendix B. Listings 
of the new and modified subroutines are given in Appendix V. Appendix D is contained 
in Volume II. 








2.0 


DISCUSSION OF METHODS 


The analytical methods utilized in the subroutines which were added to 
the SINDA routine are described in this Section. Section 2.1 describes the methods 
used for calculating temperatures of flowing fluids and Section 2.2 describes 
methods used in the pressure/flow analysis of flow networks. 


2.1 Thermal Analysis Methods 

The Flow-Hybrid method for obtaining temperature solutions was formulated 
for use with several SINDA temperature solution routines including CNFRWD, CNFAST, 
CNBACK, CNFWBK, CINDSS, HYBRID, SNFRDL, SNFRWD , CINDSL, FWDBCK, SNDSNR and STDSTL. 

The formulation included utilization of the fluid flow analysis data for the thermal 
analysis thus minimizing input and data storage requirements. The Flow-Hybrid method 
is described separately below for explicit methods (CNFRWD, SNFRDL, SNFRWD, and 
CNFAST), implicit methods (CNBACK, FWDBCK and CNFWBK), general hybrid methods (HYBRID) 
and steady state (CINDSS, SNDSNR, STDSTL and CINDSL). 

2.1.1 Flow-Hybrid Solution for Explicit Problems 

The fluid nodes temperatures are solved using the “Flow-Hybrid" solution 
method in the explicit SINDA user subroutines, CNFRWD, SNFRDL, SNFRWD and CNFAST. 

This method requires that the finite difference equations be written in the implicit 
form for the fluid lumps, while remainder of the lumps in the problem are solved 
using the explicit methods. 

The finite difference equations for the Flow-Hybrid method are as follows: 
For the fluid lump 
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where: 


For the tube lump 
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Tf = the fluid lump temperature 

Tf = the tube lump temperature 

At = time increment 
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W f = weight of fluid lump 

= capacitance of the fluid 

w = flowrate in the tube which contains the fluid lump 

•a* 

Cp - mean specific heat for the flowing fluid between the 
upstream lump and the fluid lump 



= the enthalpy of the fluid lump at temperature 

h y = the enthalpy of the fluid lump at temperature T u 

T u = the temperature of upstream lump 

HA - the convection coefficient times area 

G. t j = the conductance value from tube lump t to lump j 

w fc = weight of tube lump t 

c t = specific heat of tube lump t 
Qf = the heat absorbed by fluid lump f 

Q t = the heat absorbed by tube lump t 

If the fluid lump is the first in the tube, h u is determined as follows: 



where h ok 
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is the enthalpy of the fluid leaving tube k and entering 
fluid lump f 

is the flow rate of tube k 
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The value for T u for the first fluid lump in a tube is obtained by revere interpolation 
of the enthalpy curve at h u . The primed temperatures in equations (1) and (2) represent 
temperatures at the end of the iteration; the unprimed temperature represent these at 
the iteration start. 
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The fluid hybrid solution methods are derived as follows: 
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Solve equation (2) for T t 
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where T^-'is the intermediate tube temperature that would be obtained with no 
connection to the fluid lump. 

If equation (3) is substituted into equation (1) and simplified we get: 
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The value of in equation (3) is given by 
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Examination of equation (4) reveals two primed temperatures: T. and 

l u 

T.. . Thus, we must calculate these values prior to evaluation of equation (4). 

Z\ i 

The value of T u can be obtained if the order of calculations start with the first 
lump in the system and progresses around the systan in order, one lump at a time. 
Since the value of T^. given by equation (5) contains no primed values, its value 
may be evaluated first. Thus, the order of calculations are: 

(1) Calculate the value of all T^ using the normal explicit 
temperature calculations assuming no fluid lump convection 
exist. This is given by Equation (5). 

t 

(2) Calculate the value of all T^ in order of their position 
in the tubes starting with the first lump in the first tube 
and progressing around the system. 

(3) Update the tube temperature using equation {3) to obtain T^ . 

Of course the coefficients in equations (3), (4) and (5) are evaluated prior to 
evaluation of the equations. Methods used in determining coefficient values are 
discussed in Section 2.1.5. 


2.1.2 Fluid Temperature Solution for Implicit Problems 

The implicit user subroutines, CNBACK, FWDBCK and CNFWBK, were modified 
so that the fluid temperatures are calculated simultaneously with the other tempera- 
tures of the problem. For CNBACK, FWDBCK and CNFWBK, the fluid lump temperatures 
are calculated using the relation 
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Where all the variables are as defined for equations (1) and (2) and T t is the last 
calculated value of the tube lump temperature. 

The tube tonperatures are calculated using the normal eauations but are 
modified to use the HAT^ and HA terms as follows: For CNBACK: 
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For CNFWBK 
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The order of calculations for the implicit routines are: 

(1) Calculate the value of all using equation (6). Fluid 
flow data is utilized to obtain coefficients in the 
equation. 

(2} During the calculation in (1) the HA values and the fluid 
lump number for each tube lump are stored in the X array. 

(captured dynamic storage) 

(3) The temperatures for the remaining lumps are calculated 
using the normal calculations, except tube lump tempera- 
tures equations are modified to include the HA and HAT 
terms as shown in equations (7) and (8). 

2.1.3 Fluid Temperature Solution for General Hybrid Problems 

The HYBRID user subroutine was modified to permit calculation of fluid 
lump temperatures during the normal temperature calculations. Explicit and implicit 
lumps are determined by calculating the CSG value for each lump and comparing it with 
the input time increment. Those lumps with CSG values larger than the input time 
increment are explicit and the remaining lumps are implicit. If the tube lumps are 
all explicit, the fluid and tube lump temperatures are calculated using equations 
(3) and (4). If any of the tube lumps are implicit, the fluid lumDS are calculated 
using equation (6) and tube lumps are calculated using the following equation: 
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Where 




tt 




= the point in the iteration for evaluating heat flux 
= the conductance value from tube lump t to lump ji, 
where lump ji is implicit 

= the conductance value from tube lump t to lunp je, 
where lump je is explicit 

= the heat absorbed by tube lump t evaluated at t + ciAt 


The following calculation procedure is used: 

(1) Determine the explic ; t and implicit lumps by comparing the CSG 
value for each lump with the input computation interval. 

(2) Determine which of the explicit lumps determined in step (1) 
are interface lumps. Interface lumps are subsequently treated 
as implicit lumps. Interface lumps are any explicit lumps which 
are connected to an implicit lump. If any tube lump is treated 

as an implicit lump, all tube lumps are treated as implicit lumps. 

(3) Calculate explicit lump temperatures using the same Z G (T,-T. ) 
for all interface nodes and for adjacent lumps if the adjacent lumps 
are implicit. 

(4) Calculate fluid and tube lump temperatures for tube lumps which are 
explicit. 

(5) Calculate implicit lump temperatures. 

2.1.4 Fluid Temperature Solution for Steady State Problems 

The steady state solution subroutines, CINDSS, CINDSL, STDSTL were modified 
to include a calculation of fluid lump temperatures. The following relations are used 
for calculating the fluid and tube lump temperatures: 
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The following order of calculation is used: 

(!) Calculate the fluid lump temperatures. 

(2) Calculate temperatures for all other lumps using the 
normal SINDA temperature equation, except tube lump 
temperature equations are modified to include the HA 
and HAT terms. 

2.1.5 Coefficient To Temperature Equations 

A brief description of the methods used to obtain the convection heat 
transfer coefficient, H, and the mean specific heat for a flowing fluid Cp (Equation 1) 
is considered to be of value. They are discussed separately below. 

C onvection Coefficient Determination 

Several methods are available to the SINFLO user for determine the heat 
transfer coefficient, H. The different options are available for each fluid lump and 
are specified by supplying values for FI, the eighth value of the type array (see 
Section 4.2). When FI is real, the programmed equations for flow in a tube are used 
to obtain H. Using this method, the flow regime is assumed to be laminar when the 
Reynolds number is 2000 or less. For this regime the convection heat transfer coef- 
ficient is calculated by: 


h = k 
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where : k 

- 

thermal conductivity 

D 

= 

hydraulic diameter to flow 

X 

= 

distance from tube entrance 

Re 

- 

Reynolds number 



• 

4 m 



FP 

m 

= 

flow rate of fluid 

f* 

= 

viscosity of fluid 

P 

= 

wetted perimeter of fluid flow passage 

FI 


An input factor for modifying fully 
developed flow 

F2 

= 

An input factor for modifying developing 
flow 


Equation (12) is a curve fit obtained by VSD to approximate the Graetz solution to 

V 1 

flow in a tube for values of greater than 0.001. The convection heat transfer 

coefficient for flow in a tube in the transition flow regime (2000 < Re < 6400) is 
approximated by the following relation: 

h = | [0.116 (Re 2/3 - 125) (Pr) 1/3 ] (13) 

This relation was derived by Hausen and holds only for fully developed flow. The 
relation used to determine h for turbulent flow (Re 6400) is the following 

h = .023 | (Re)* 8 (Pr) 1/3 (14) 

If FI is the integer 1, a more general option is used for determining the 
convection heat transfer coefficient. A curve of St(Pr) 2 ^ 3 vs Reynolds No. is inter- 
polated to obtain the value of St(Pr) 2 / 3 . That is, 

St(PR) 2/3 = F(Re) (15) 

where: St = Stanton number 

Nu 

Re Pr 
h 

CpV 
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V = Average fluid velocity 

F (Re) = An arbitrary function of Reynolds number which the 

user inputs as a table (identified by F2) 

The heat transfer coefficient is calculated by 

h = | Me) Re(Pr) 1/3 (l6> 

If FI is input as the integer 2 the convection heat transfer coefficient is obtained 
by direct interpolation of a curve of heat transfer coefficient vs flowrate which is 
identified by F2. 

Mean Specific Heat 

The method for obtaining the convective term (wCp) in equation (1) was 
modified with SINFLO to get better accuracy. Rather than just using the lump specific 
heat, a mean value of specific heat was obtained as follows. 

The mean specific heat for a fluid going from upstream temperature T^ 
to fluid lump temperature T f may be obtained by integrating: 




(17) 


Where h is the fluid enthalpy 

Equation (17) is used to obtain the mean specific heat for equation (1). The enthalpy 
curve is required as input. 
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For temperature lumps that have multiple upstream lumps, such as mixing 
junctions, the value of h u is determined as follows: 



(18) 


and the upstream temperature, T^, obtained by reverse interpolation of the enthalpy 
curve. Equation (17) is then applied using the values of h u and thus obtained. 

2.1.6 Heat Exchanger Analysis 

Five subroutines have been written to facilitate the thermal analysis 
of systems containing heat exchangers. These are HXCNT for analysis of counter 
flow heat exchangers, HXPAR for parallel flow heat exchangers, KXCROS for cross 
flow heat exchangers, HXEFF for any heat exchanger with an input effectiveness, 
and HXCOND for condensing heat exchangers. These subroutines calculate the outlet 
temperatures of two sides based upon the inlet temperatures and heat exchanger 
effectiveness. The relations used for calculating effectiveness are described below. 


2. 1.6.1 Counterflow Heat Exchanger 

Subroutine HXCNT calculates the heat exchanger effectiveness using the 
relation from Reference 1 for counterflow heat exchangers. That is. 



Where 


effectiveness 


( 19 * 


LfA = overall effectiveness 

(MC) S = mass, specific heat product for the side with 
the smallest MC 

(MC)j = mass, specific heat product for the side with 
the largest MC 
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The limiting cases for this relation are: 


(!) When (MC)_/(MC) = 0, 

b 1 

e = 1-e -UA/(MC) s 

(2) When (MC) S /(MC) I = 1 

UA 

TO? _ UA 

e " 1 + UA (MC)Y + UA 

TOY 


Using the effectiveness as calculated by the above method, the outlet tempe**- mres 
are calculated as follows: 

1. For the side with the smallest MC, (MC)s : 


Tout s = Tin s - e (Tin s - Tin^) 


( 20 ) 


2. The enthalpy of outlet for the side with the large MC is 
then calculated by 


hout^ 


hin^ + ( hi n $ - hout s ) 


w. 


w. 


( 21 ) 


where: 


hout^ 

hin £ 
hi n s 

hout s 


enthalpy of the outlet for the side with the large MC 
enthalpy of the inlet for the side with the large MC 
enthalpy of the inlet for the side with the small MC 
enthalpy of the outlet for the side with the small MC 
flow rate of the side with the small MC 
flow rate of the side with the large MC 


3. Tout^ is obtained by reverse interpolation of the enthalpy curve 

at hout^. 


2.1. 6. 2 Parallel Flow Heat Exchanger 

Subroutine HXPAR calculates the heat exchanger effectiveness using the 
relation for parallel flow heat exchangers^ which is: 
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The limiting cases are 

(1) When (MC) S /(MC) / - 0, 

« . ! _ e “ UA/(MC) s 

(2) When (MC) s /(MC)i = 1.. 

, UA 
' 2 THC) S 
<■ = 1 - e 


( 22 ) 


2.0 

The heat exchanger outlet temperatures are then calculated using the 
method described for HXCNT. 

2. 1.6. 3 Cross Flow Heat Exchanger 

Subroutine HXCROS calculates the effectiveness for cross flow heat 
exchangers using one of the four relations below depending upon mixing of the streams. 
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The heat exchanger outlet temperatures are calculated using the method described 
for HXCNT. 

2. 1.6. 4 User Supplied Effectiveness 

Subroutine HXEFF was written to perform heat exchanger thermal analysis 
with a user supplied effectiveness. The effectiveness may either be supplied as 
a constant or as an array number which gives the effectiveness as a bivariant 
function of the flowrates on the two sides. The outlet temperatures are then 
calculated using the method described for HXCNT. 

2. 1.6. 5 Condensing Heat Exchanger 

Subroutine HXCOND was written to analyze a condensing heat exchanger. 
The effectiveness may either be supplied as a constant or as a trivariant 
function of humidity, flow rate of the gas, and flow rate of the coolant. The 
outlet temperatures are calculated as follows: 



TG out 

' TG in * s < TG in ‘ TC in> 


where: 

TG out 

= temperature of the gas 

out of the heat exchanger 


e 

- effectiveness 



TG 1n 

= temperature of the gas 

into the heat exchanger 


Tc in 

= temperature of the coolant into the heat exchanger 


1 , r j ■ 
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The saturation pressure is given by 


TG .-500 
(19.3 — 2^ ) 


PB out ” - 1217e 


out 


where: PB 0U t = saturation pressure of the gas 

And the outlet humidity is 


XMIM0 • PBOUT 

^ = 

P - PBOUT 

Where \p = humidity Ojj- n > > 0) 

XMIM0 = molecular weight ratio 
P = total gas pressure 

The flow rate of the liquid is 


where: 


*1 - *g - *out J 

&£ = flow rate of the liquid 

ft = flow rate of the gas 

9 


The enthalpy of 



the coolant out of the heat exchanger is 

[ «hg 1n + hg out ) * Q > + ■ XLAM ] 

= hc in + 7T^ 


where: XLAM 


latent heat of vaporization 


The outlet temperature of the coolant is obtained by reverse interpolation of the 
enthalpy curve at hc QUt . 
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2.1.7 Cabin Analysis 

A subroutine has been written for use with SINDA which will give 
the user the ability to perform thermal analyses on cabin air systems including 
condensation on the walls and a vapor mass balance. The cabin heat transfer 
and condensation analysis involves the two-component flow of a condensible 
vapor and a non-condensable gas, with condensation of the vapor occurring 
on surfaces in contact with the fluid. Two problems of this nature have 
been studied extensively. 

1. Condensation on, or evaporation from, a surface over which 
a free stream of fluid is passing. In this case, for rela- 
tively low mass transfer rates, the fluid properties can be 
assumed to be constant. 

2. Dehumidification of a confined fluid stream by a bank of 
tubes. In this case there is a marked change in the 
temperature and vapor content of the fluid, and the de- 
tailed deposition of the condensate is not of primary 
interest. This type of analysis is usually handled or: an 
overall basis similar to heat exchanges effectiveness 
calculations. 

The following additional assumptions have been made with respect 
to the cabin atmospheric conditions. 

1. The heat of circulation in the cabin is sufficiently high 
that the temperature and humidity are effectively the same 
throughout the cabin. 

2. The velocity at all points where heat transfer and/or con- 
densation can occur is known, and is proportional to the 
total mass flow rate in the cabin. 

These assumptions' make it possible to calculate the heat and vapor 
balance in the cabin for the entire volume as a unit, and to solve the heat 
transfer and condensation equations at each node independently of the other 
nodes. 



Cabin humidity can be determined from an overall vapor balance in 
the cabin. The total vapor in the cabin at the end of an iteration is: 


W V =W V + W v in - W v out - SW L 

Where W v = mass of vapor in cabin at end of iteration i 

W v 1_ = mass of vapor in cabin at start of iteration i-1 

W v j n = mass of vapor flowing into cabin during iteration i 

W v 0U |. = mass of vapor flowing out of cabin during iteration i 

2 W|^ = mass of vapor condensed during iteration -1 


W v -j n is determined from the known conditions of the gas flowing into the 
cabin. 


Where 


in = iti in 

_ 


^in 

1.+ iPin 


m in = mass flow rate into cabin 
& in - specific humidity of gas flowing into cabin 
= time increment 


It is assumed that an equal volume of gas is flowing out of the cabin. 
Then » 

w v out = m out f , V I 

L 1 + <A C J 

Where & = specific humidity in the cabin (at the end of the 

previous iteration) 

and rh out = in in [ft /Pin] 

Where p„ = cabin density 

c 

pin = density of gas flowing into cabin 


The condensation term 2Wj_ is determined from the calculations for 
the individual nodes as described below. The properties of the cabin 
atmosphere are determined from the calculated value of Wy. The vapor pressure 
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1 


I 


in ilia cabin is 


Pv 



R v T c 


Where 


Vc 

Rv 

T c 


cabin volume 
gas constant 
temperature of cabin gas 
vapor pressure 


Assuming that the cabin pressure 
partial pressure P a is: 


= Pc * P v 


is a constant, the gas 


and W a = 

R a T c 

Where W a - mass of non-condensible gas in the cabin. 

Now the new value of specific humidity in the cabin can be 
determined by 



Where 


The properties of the atmosphere can now be determined by 

Me = XHq + tyc M y 
X + ^C 

Cpc = Cpg -f^cCpv 
1 + c 

k r = Xkq + <Ack v 

— X + WI 

Pc = w v + W s 
’ h 

M = viscosity 
Cp - specific heat 

k = thermal conductivity 
X ~ molecular weight ratio, M v 
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and all values are evaluated at L, . Cabin temperature T c can be 
determined by a heat balance on the cabin atmosphere. 

Tc = Tc ’ 1 + 111 1n C P c (Tin - Tc 1 ' 1 ) - £Qi 

(Wy + W/\) Cpc 
i -1 

Where T c = T c after previous iteration 

Ti n = temperature of gas flowing into cabin 
2Q L - net heat loss to cabin lumps 

The heat transfer between the cabin atmosphere and the tube 
and structure lumps in the cabin is defined by: 

QL: = hA Li [ Tc - Tj_i ]Ar 


Where 

h 

= heat transfer coefficient 
= heat transfer area of lump 



Tli 

= temperature of tube lump 




= time increment 


Using the Colburn-Chi Iton heat transfer-mass 
the condensation (or evaporation) at the tube lump is 

; transfer analogy, 
determined by: 



= ^m A Li [ p v “ p wi] A r 


Where 

W Li 

= condensation on wall, lb. 



^lil 

= mass transfer coefficient 



p wi 

= vapor pressure at T|_j 


is 

The latent heat addition to the lump due to 

this condensation 


h 

<3 

AW U X 


Where 

A = 

latent heat of vaporization 


The vapor pressure P w i can be determined by a relationship 
derived from the Clausius-Clapeyron equation and the perfect gas law 
(Appendix K of Reference 3). 
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Where P 0 is known vapor pressure at a reference temperature 


Three methods are available for determining mass and heat transfer 
coefficient. For tube lumps the equations from Reference 1 for gas flowing 
normal to the tube axis was assumed. Three different equations are used 
depending on the value of the Reynold's number. 

Nu * 0.43 + .533 (Re)* 5 (Pr)* 31 Re < 4000 

Nu = 0.43 + .193 (Re)* 618 (Pr)* 31 4000 < Re < 40000 

Nu = 0.43 + .0265 (Re)* 805 (Pr)* 31 40000 < Re < 400000 

These equations were derived for an air-vapor mixture, but should 
be relatively accurate for other similar gases. The Nusselt and Reynold's 
numbers in the equations are defined using the tube diameter for the 
characteristic dimension, and the velocity in the Reynold's number is input 
at each lump and ratioed to the total cabin atmosphere flow rate. 

Wc 

v\ - yio 

Wco 

• 

Where Wco = nominal cabin atmosphere circulation rate 

^io = velocity at lump at Wco 

W c = circulation rate at time of calculation 

The second option assumes flat plate flow for cabin wall lumps. 

In this case the heat transfer coefficient, for laminar flow, varies along 
the plate. Hence, direction of gas flow and the location of an assumed 
leading edge must be assumed. The equation for flat plates from Reference 1 is: 

N u = 0.332 Re* 5 Pr 1/3 

where the Nusselt and Reynold's numbers are local values and are 
defined by the distance X from the assumed leading edge. For a wall lump of 
length Li which is located a distance L-j 0 from the assumed leading edge, the 
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average Nusselt number can be defined as: 

N u = 0.664 Pr 1/3 j^Re^* 5 - (Re 0 )‘ 5 

Where Nu is defined by L] 

Re 0 is defined by L] 0 
Re] is defined by Li 0 + L-j 


The third option is a direct user input for convective heat transfer coefficient. 

For the determination of mass transfer coefficients, the same 
equations which were used for heat transfer coefficient can be used with the 
Sherwood number substituted for Nusselt number and Schmidt number for Prandtl 
number. However, if the diffusion coefficient for the cabin is approximately 
equal to thermal diffusivity, the Sherwood number is equal to the Nusselt 
number and the mass transfer coefficient can be determined directly from the 
heat transfer coefficient. That is: 


Sh = Nu 

KmRTgX = h x 

T iT 

If D =* a then 

Km = hD 

apCpRTg 

Km = -JL... 

Cp P c 


(28) 


Equation (28) is the Lewis relationship (Reference 1). For a mixture 
of oxygen and water vapor characteristic values are .866 for the diffusion 
coefficient, D, and .879 for thermal diffusivity, a, so the relationship 
should be valid. 

For cabin tube and wall lumps the values for AQ^ and AQx] are 
added to the basic heat balance equation for these lumps. Values for AQu 
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are summed for all participating lumps for input to the cabin atmosphere heat 
balance. Values for AWLj are also summed for all lumps for cabin humidity 
balance, and the value for total water condensed on each lump WL.. is main- 
tained. 

If the rate of evaporation or condensation is high it would be 
possible for the cabin humidity to change significantly during a single 
iteration. This could lead, for example, to overestimating condensation by 
assuming that the humidity is constant in the calculation. A test of the 
approximate vapor pressure in the cabin at the end of the iteration is made, 
and the condensation or evaporation at any lump is reduced, if the sign of 
the AWL^ term is changed. A value Wy' is calculated by: 


V = W v i_1 - SW L i 


P v l = W v ' 


T4TVT 


R v Tg 


Then for each lump if 


p v 1 “ p wi < o 

p v * p wi 


a new value of AW^i is calculated by; 


A Wjj - A 


r 

«u r 


The new values of AW^ are now again summed for the new value of 
ZAWl for establishing cabin humidity for the next iteration. A test is also 
made to assure that Ww' is never less than zero. 


2.2 


Fluid Flow Analysis 

Subroutine FL'? C QL was written as a SINDA user subroutine to provide 
the ability to perform fluid pressure/flow analysis for flow of an incom- 
pressible fluid in tubes. The fluid flow analysis of FLOSOL is integrated 
with the thermal analysis capability so that the temperature dependence 
of properties is included in the pressure balances. FLOSOL is called from 
the VARIBLES 2 user logic block. 

FLOSOL performs a pressure-flow balance on a general flow network 
including the following effects: 

(1) Friction pressure drop 

(2) Orifices and fitting type pressure losses 

(3) Valves 

(4) Pumps 

(5) Incoming flow sources at any pressure point in the system 
The user describes the flow model to the subroutine by supplying the tube 
network connections and information concerning fluid properties, flow 
geometry, temperature model lumps, orifices, valves and pumps. Using this 
information, the subroutine determines the flow distribution required to 
satisfy (1) the conservation of mass at each node point and (2) equal 
pressure drops across tubes in parallel. The model used to describe the 
flow system and the analytical methods for determining the solution are 
described below. 

2.2.1 Overall Flow Model Description 

A flow problem may be analyzed with FLOSOL, simultaneously with a 
thermal analysis, so that the flow solution is continually updated based 
on the thermal conditions. To perform a flow analysis, the user must in- 
put a mathematical model of the flow system. The flow system is assumed 
to consist of a set of interconnected tubes such as the example shown in 
Figure 1 , which consists of two radiator panels, each containing four 
tubes and connected so that they flow in parallel. 

For clarity the following definitions are made at this point: 
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G> Tube Numbers 

XX Pressure Nodes 


FIGURE 1 FLOW SYSTEM SCHEMATIC 




(1) A tube Is any single length of pipe between two pressure 
nodes. A tube "contains" fluid temperature nodes and may 
contain as many of these as required. 

(2) A pressure node is located at each end of a tube. As many 
tubes as desired may be connected at a node junction and 

a node must exist at the junctions of two flow pipes. 

We must make a mathematical model to describe the fluid flow information to 
the computer. The information required consists of: 

(1) Identification of the pressure node numbers 

(2) Identification of the tube numbers and the two pressure 
nodes connected by tube 

(3) The fluid temperature nodes contained in each tube 

(4) The flow geometry for each temperature fluid nodes 

(5) The number of "head losses" for items such as orifices 

(6) Fluid property information 

(7) Valve connections and characteristics 

(8) Pump characteristics 

To build a flow mathematical model, a schematic of the flow system is needed. 

As shown in Figure 1, the pressure nodes and tubes may be superimposed on 
the schematic. It is also helpful to impose the fluid temperature lump numbers 
for each tube. 

To facilitate speedy analysis on a general flow problem, provisions 
have been made for the user to divide the flow system network into subnetwork 
elements. For example, the flow system shown in Figure 1 could be divided 
as shown in Figure 2 . Tubes 23 and 24 are added in the main network as shown 
in 2(a) to replace subnetwork elements 1 and 2. The subnetwork elements 1 
and 2 which are shown in Figures 2(b) and 2(c) are then input as separate 
network elements. This type of subdivision allows the solution to be obtained 
by solving two sets of 6 simultaneous equations and one set of 8 equations 
rather than the original set of 16 simultaneous equations. This type of sub- 
division has been found to enhance the solution speed and accuracy for pro- 
blems with a large number of nodes. 

In summary, the pressure/flow solution is obtained by the following 

sequence: 
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(1) The flow resistance is obtained for each fluid temperature 
lump in each tube including the effects of friction, orifices, 
and fitting type losses. 

(2) The flow conductor value is obtained for each tube by sunming 
all the resistances of the fluid lumps in the tube, adding 
the valve and user supplied resistance to the sum, and in- 
verting the resistance. 

(3) A set of simultaneous equations is set-up and solved for each 
main system and subnetwork to obtain the pressures. 

(4) The flow rates are then calculated. 

A detail discussion of each element in the above sequence is described in the 
following subsections. 


2.2.2 Tube Conductor Determination 

The value of the flow conductor is determined for each tube by 
first calculating the flow resistance for each temperature fluid lump con- 
tained in the tube, summing these resistances up to obtain the flow re- 
sistance of the tube and inverting the tube resistance to get the conductance. 
Flow conductance is defined by the relationship 



= SF ij' p i - p j) 


(29) 


Where &jj = flow rate between pressure nodes i and j 

GFjj = flow conductance between nodes i and j 

Pi = pressure at pressure node i 

g = pressure at pressure node j 

The flow resistance for each lump is then 


R 


. 1 

K " GT ~ 


ft 


Where R^ = flow resistance for lump k 
APj, = pressure drop for lump k 


But AP^ is given by 
AP k 



\ w 2 

) 


(30) 




28 



Where = the friction factor for lump k 
ffc = the friction factor coefficient 
= the lump length for lump k 
D = the lump hydraulic diameter for lump k 
K = the dynamic head losses for lump k 

ft = the flow rate 
g c = the gravitational constant 
= the fluid density for lump k 
A = the flow area 
The flow resistance is then given by 

R k = ( f k ffc + 

Two options are available for obtaining the friction factor, fy. These are 
(1) internal calculations for all flow regimes and (2) internal calculation 
for laminar flow and obtained from a table of f vs Re (where Re is the 
Reynold's number) for transition and turbulent flow. For the first option 
the internal calculations for the three flow regimes are: 

Laminar Regime : Re^ < 2000. 


K 


ft 


2g c p k A 


(31) 



Where f^ = friction factor for lump k 
Rej^ = Reynolds number tor lump k 


(32) 


Transition Regime : 2000 < Re^ < 4000 


.2086082052 - .1868265324 
+ .06236703785 


Re. 


1000 


(33) 


Re k I 2 * .0065545818 


1000 


n * 

i 


Re k 

Tooo“ 


2 ? 


Turbulent Regime : Re^OOO 



Equation (33) for the transition regime is a curve fit between the 
laminar and turbulent regimes which was derived to match the two curves in 
a continuous manner. It is merely an arbitrary curve in this undefined 
region. A curve of the friction factor vs Reynold's number given by the 
above relations is shown in Figure 3. 

The second option for friction factor uses equation (32) for the 
laminar regime and a user input curve of f^ ve Re for the other regimes. 

The options available for input of the dynamic head loss, (»( , include 
(1) an input constant or (2) a tabulated curve of ft vs Re. 

To obtain the conductance for each tube, the flow resistances for 
all the lumps in the tube are added and then inverted, giving 

GF,. =J 

j Z R k < 35 > 

k K 


2.2.3 Valve Analysis 

Two methods have been included in the FLOSOL subroutine for modeling 
valves. These are (1) pressure drop through the valve is included in the system 
flow balance and (2) valve position is used as a fraction for splitting the flow. 
The first method uses the following equation to characterize the pressure drop 
through each side of the yalve: 


AP 


t[|3 2 


where AP 
E 
W 
X 


valve pressure drop 

valve pressure drop factor (user input) 

flow rate through the side of the valve under consideration 

the fraction of the valve opening 


The second method for modeling valves uses the valve position as a fraction 
for splitting the flow into the valve according to the following equation: 


w. 


* X i*in 
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Re x 10" 3 (curve 1) 

• Re x 10' 4 (curve 2) 

I FIGURE 3 FRICTION FACTOR VS REYNOLDS NUMBER 
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where W. = flow rate out side 1 of the valve 
X.. - valve position of side i 

*in = ^ ow ra ^ e the va ^ ve 

The pressure drop through the valve is not included in the system flow balance. 

The valve pressure drop factor controls which method is used. The first 
method is used when the valve pressure drop factor is greater than 0. To specify 
the second method the user must input a 0 for the valve pressure drop factor. For 
either method the value of X must be greater than 0 and less than 1. 

Three basic types of valves are available in FL0S0L for either the pressure 
drop or the flow splitting formulations which give different characteristics for the 
dynamics of the valve position X. These types are: (1) Rate Limited; (2) Polynomial, 
and (3) Shut-off. 

A number of variations are available for each valve type. For instance, each 
of the above may be either one sided or two sided. If a valve is two sided, the 
valve position of side 2, X 2 , is related to that of side one by 

x 2 = 1.0 - X 1 

If the valve is one sided, either side one or side two may be used. Provisions 
are included for a valve time constant to be included with the polynomial valve. 

The methods used to obtain the valve positions for each of the three 
types are described below. 

2.2. 3.1 Rate Limited Valve 

The valve position for the rate limited valve is obtained by an 
approximate integration of the valve rate of movement, X. X depends on 
the temperature difference between the valve control set point temperature 
and the sensor temperature as shown in Figure 4. With this characteristic, 
the valve has no movement as long as the valve temperature error, ATJs 
within the dead band. Outside the dead band, the velocity of the valve in-* 
creases linerarly as the error increases to a maximum rate, X max. The dead 
band, rate of velocity increase, dX/d{AT), and the maximum velocity are 
control,! ed by user input. 

The relations used to obtain the valve positions are as follows: 


X 1+1 - X 1 + (* ,+1 ) (Ai- ) 


(36) 


Where X i+1 
X 1 

X 1 * 1 

Ar 


valve position at iteration i+1 
valve position at iteration i 
valve velocity at iteration i+1 
the problem time increment 
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Valve Velocity, X = Fraction Open per Unit Time 
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FIGURE 4 RATE LIMITED VALVE OPERATION 


The valve position is limited by 

X min ^ X i+1 < X max 

Where X min and X max are input limits on the valve position. 
The valve velocity, X x *^ , in equation is given by: 



i 1+1 

= 0 

if 

j Tsen 

- Tset| £ 

where: 

x i+1 

• 

= [Tsen-Tset-Tdb] 

if 


> Tset + Tdb 


i 1+1 

* 

= [Tsen-Tset+Tdb] 

if 

Tsen 

< Tset - Tdb 


Tsen = sensor lump temperature 

Tset = set point temperature 

Tdb = valve dead band temperature 


The valve velocity is limited by 

* • {4.I • 

Xmin ^ X < Xmax 


After the valve position for side 1 is obtained from equation the side 

2 position is obtained from x*> = 1*0 - X] 

2. 2. 3. 2 Polynomial Valve 

The polynomial valve determines the steady state valve position 
as a 4th degree polynomial function of the temperature error between the 
sensor lump and the set point. A valve time constant is then applied to 
determine how far between the previous position and the new steady state 
position the valve will move. The steady state position, X ss , is given by 

X ss = A o * A i AT + A 2 Aj2 + A 3 Aj3 + A 4 at4 

Where AT = Tsen - Tset 

Tsen = the sensor lump temperature 
Tset = the set point temperature 
A o s A 1 * A 2* A 3* A 4 = ln P ut constants 
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The valve position, X i+1 is then determined by 


,i+l 


= X ss + < x1 - x ss> e ' 4r/rc 


Where X 1+1 
X 1 


r c 


valve positon at iteration i+1 
valve position at iteration i 
problem time increment 
valve time constant 


The valve position for side 2 is given by 

x 2 = 1.0 - x 1 

where is given by equation (37) 


(371 


If one desires to eliminate the effect of the time constant (and thus, give 
the valve an instantaneous response), a value for x c shouid be input which 
is small compared to the time increment, At. Also, either a constant value or 
a temperature lump number may be specified for the set point to permit use 
of the valve for proportioning between two sides. 

2. 2. 3. 3 Shut-off Valve 

For side 1 of a shut-off valve the valve position decreases from 
X max to Xmin when the temperature of the sensor lump drops below the specified 
"off" temperature* T 0 ff, and increased from X m -| n to X max when the sensor lump 
exceeds a second specified temperature, T 0 n* Ton must be greater than Toff- 
Side 2 works in reverse of side 1. The valve position increased from Xjy,-} n to 
Xjnax when the sensor temperature drops below the specified T on and decreases 
from X^x to X m -j n when the sensor lump increases above the off temperature, 
T 0 ff. For side 2, T 0 ff must be greater than T on . Note that, if the shut-off 
valve is a two sided valve with both sides active, the valve is a switching 
valve. 

2. 2. 3. 4 Valve Flow Resistance Calculations 

The valve pressure drop on side one is assumed to be given by: 

AP = E [|] 2 (38) 




f 
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Since flow resistance is AP/w, the valve flow resistance is given by 


R. 


Eft 

7 


( 39 ) 


This value of flow resistance is calculated and added to the other flow resistances 
of the tube prior to performing the operation in equation (35) to find the value 
of the flow conductor for the tube. 

Valves may be either one way or two way, i.e. , be one tube or two tubes 
at the outlet. If only one tube exists on the valve outlet the flow resistance is 
calculated using equation (39) above. If a second tube exists, the resistance 
on side 2 is given by 



E 2*2 

O-x ) 2 


( 40 ) 


2.2.4 Pressure-Flow Network Solution 

As previously stated, the user may subdivide a system flow net- 
work into a main network and subnetwork elements. The elements which are 
subnetworks to the main network may also contain subnetwork elements but 
the subdivision can go no lower than two levels. 

After the flow conductor values havq been obtained by the methods 
described in Sections 2.2.2 and 2.2.3 a set of simultaneous equations are set 
up and solved for the main system and for each subnetwork. The subnetwork 
elements are all solved first and then, their equivalent flow conductor value 
is calculated. The value is inserted in the main system network and the 
system solution is obtained. The procedure is repeated until the problem 
is balanced. 

A set of simultaneous equations are obtained by conservation of mass 
at each pressure node for each network and subnetwork. For any node i the 
conservation equation can be written as follows: 
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E* 


out 


V w. = 
^ in 


= 0 


Let 


w _ w 
in" i 

nc 


and Iw out = ^ SF jj 1 Pj - P, 1 


0 = 1 


Then equation (41 ) becomes 


J=1 


(41) 


E GF ij !P j - p i ! - *i ' 


* 0 


1-1 .n (42) 


Where GF.r • = flow conductor between pressure nodes i and j 

I J- 

= pressure at node i 
- pressure at node j 
= flow rate added at node i 
= number of pressure nodes in the subnetwork. 
Tite above equation is a set of n simultaneous equations for P array. 
Pressure in the system or subsystem may be set at a specified level but 
the last (outlet) node must be specified. Equation (42) may be written 
in matrix form as: 


GF « 

p i 

P .J 

W. 

1 

n 


Where 


GP = C 


( 43 ) 


£ GF ij " GF 12 “ GF 13' 
’ GF 21 E GF 2j " GF 23 


n-1 


“l + GF 1n p n 
h + GF 2n P n 


w , , GF , -i P 
n-i + n-1 1 n 


P n is the specified pressure. The above equations are solved for pressures 
at each point in the system and flow rates are then calculated for each 
tube by: 




GF u t p i 


p j> 


Since the coefficient matrix given by equation &3 ) is symetric and positive 
definite the efficient square root or Symmetric Cholesky method was pro- 
gramed to obtain the solution. This method is more accurate and faster 
than any other methods studied for this application. 

Since the flow conductors are functions of the flow rate, the set 
of equations given by (43) are solved numerous times on each temperature 
iteration with a new set of GFjj values for each solution. The iteration 
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process continues until the change in the flowrates is within some user 
specified tolerance before proceeding to the next iteration. 

2.2.5 Pump and System Pressure - Flow Matching 

Concurrent with iterating the system flow equation to solution 
on each temperature iteration the overall system pressure drop and flowrate must be 
matched to a pump characteristic. Several types of pump characteristics 
are available to the user as options. These are (1) system flow rate 
specified as a constant, (2) system flowrate specified as a known function 
of time, (3) pressure drop specified as a function of the flowrate in a 
tabulated form and (4) pressure drop specified as a function of flowrate 
with a fourth degree polynomial curve. 

The first two options require no balancing of the pump with the 
system. Balancing is required for options (3) and (4) and iterative 
procedures have been devised to obtain the solution of the pump curve 
to the system characteristics with as few passes as possible through the 
system pressure/flow balancing loop for these options. The procedures 
used for these options are described below. 

2. 2, 5.1 Tabulated Pump Curve Solution 

The matching of a tabulated pump pressure rise/flow characteristic 
to the system pressure drop/flow characteristic is accompli shed *by the 
following procedure. See Figure 5 to aid in understanding the procedure. 

Step 1 : The initial flowrate, ft], at the system inlet is 
established either from user input on the first 
iteration or the system flow of the previous iteration 
for subsequent iterations. 

Step 2 : Using ft], a solution to the flow network is obtained 
using the methods described in Sections 2.2.2, 2.2.3 
and 2.2.4, Following this solution, AP] is available 
establishing point 1 on the true system characteristic 
curve shown in Figure 5. 

Step 3 : Obtain an equation for the straight line approximation 
of the system characteristic (line 0, 1 for the first 
pass, line 1, 2 for the second pass, etc.) 

AP S = C W s + D 
oo 


Pressure Rise or Drop* AP 







5 tt>^ r iiT 


where C 



AP 


0 



D = AP 0 - 


AP 1 - AP 0 * 
*1 ■ *0 c 


AP S> & s are the system pressure drop and flowrate 
values given by the approximate equation 
AP] , w-j are the latest values for system pressure 
drop and corresponding system flowrate 
APq, tig are the values for system pressure drop and 

corresponding system flowrate for the previous 
pass (These values are zero for the first pass) 

Step 4 : Obtain the equation of the line connecting points ai 

and b] which is an approximation of the pump characteristic. 

(1) Two points are determined on the pump characteristic 
curve: 

(a) interpolate the tabulated character istic at ft] 
to obtain AP a ] (See Figure 5) to locate point ai 
at W] , AP a ]. If W-j is greater than fynax 5 set W*] 
equal to w max and AP a ] equal to zero. 

(b) reverse interpolate the tabulated characteristic 
at APl to obtain to locate point b] on the curve. 

If AP] is greater than AP max , AP] is set to AP max 
and ti/^] is set to zero. 

(2) Determine the coefficients A and B for the equation 


where A 


AP p = Artp + B 

AP 1 - Ap al 
W bl " W 1 


B 



AP 1 - AP al 
w bT W 1 


*1 






APp, *p are the pump pressure rise and flowrate as given 
by the approximation. 

Step 5 : Solve the approximate equations obtained in Steps 3 and 4 
to obtain an approximate solution to the system 
characteristic and the pump characteristic (Point N) as 
follows : 


fo - D - B 
N ~ A - C 


APj^j = Afog + B 

Step 6 : Check the tolerance below where fo n _-| is the previous 
&N (fo-| for the first time through) 


'N ■ W N-1 . 
*N- 


(1) If the above inequality equation is not satisfied repeat 
steps 4 through 6 substituting % for ft] and 

APpi for APi 

(2) If the inequality is satisfied the point SI (Figure 
5) has been located. Continue with step 7. The 
final flowrate is i/z 

Step 7 : Check the following tolerence 


Is 




< TOL* 


(1) If the above inequality equation is not satisfied, 
repeat steps 2 through 7 using the value of fog 
for fo] . 

(2) If the Inequality is satisfied, fog is the solution 
flowrate. 


* 


TOL is the input pressure solution tolerance 


r 


1 
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2. 2. 5. 2 Polynomial Pump Curve Solution 

When the user describes the pump curve with a polynomial curve 
fit, the pump characteristic is described by the relation 

AP p = A 0 + A 1 * + A 2 W 2 + A/ + A/ 

When this option is used, the procedure for matching the pump characteristic 
to the system characteristic is identical to that described 
in Section 2. 2. 5.1 for the tabulated pump characteristic except Steps 4 and 
5 are replaced with the following: 

Step 4 : Obtain the coefficients of the 4th order equation 
to be solved 
Since: 



AP = C ft + D (C and D are obtained from Step 3) 

A p p = 'b + Vp + A 2*p + A 3*p + A 4*p 
The solution occurs when 
AP S ~APp 

Then the equation for ft^ is 

(V D) + (Ai“C) ft N + A 2 ft 2 + A 3 ft^ + A 4 ftjJ = 0 

Step 5 : Solve the equation for ft^ using the Newton-Raphson 
Method of solution for a fourth order polynomial 


The remaining steps are identical to that given in Section 2. 2. 5.1. 


3.0 


SINDA ROUTINE MODIFICATIONS AND ADDITIONS 

This section describes the actual modifications and additions to the 
SINDA program. Preprocessor changes are discussed in Section 3.1. Processor 
changes are discussed in Section 3.2. The specific changes can be seen in the 
listings in Appendix C. 

3.1 Preprocessor Modifications and Additions 

The preprocessor required four interface points to include the FLOW 
DATA block. The first point was in subroutine CODERD prior to the processing of the 
CONSTANTS DATA where a call is made to subroutine FL0W1 which reads and interprets 
the FLOW DATA cards. The second point is also in subroutine CODERD immediately 
following the processing of the ARRAY DATA where subroutine FL0W2 is called to pro- 
cess the FLOW DATA. The last two interface points involve setting up arrays in 
labelled common blocks. Subroutine FLOCOM does this when called from subroutine 
GENLNK and subroutine PRESUB. 

Subroutine PSEUDO was modified to allow a node to not have any connections 
in the BCD 3CONDUCTOR DATA blcck. This was necessary since the tube lump to fluid lump 
connections are internally generated and not defined in the BCD 3C0NDUCT0R DATA block. 

3. 2 Execution Routine Modifications 

The requirements for integration of the Fluid Hybrid Temperature solution 
subroutine, FLUID, with the SINDA temperature solution subroutines were minimized 
at the interface point. A labelled corrmon block which contains a code for flow problems 
was added to each of the temperature routines. When the code is tested for a flow 
problem, subroutine FLUID is called to calculate the fluid lump temperatures. In the 
explicit routines the tube lump temperatures are also calculated. In the implicit 
routines the convection conductance and adjacent fluid lump number are stored for 
the tube lump and after returning to the mainline temperature routines, the convection 
conductance is included in the calculation of the tube lump temperature. The SINDA 
routines were modified so that the node with the minimum natural time increment would 
not be a tube lump which does not have any connections from the conductor data block. 
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4.0 FLOW DATA BLOCK INPUT FORMAT 

This section describes the input format for the SINFLO input data. The 
lumps referenced in this section must have been entered in the NODE DATA block. 

The fluid lumps must be entered as boundary nodes and the tube lumps must be 
entered as diffusion nodes. 

The SINFLO input data for the fluid systems are supplied by the new data 
block headed by "BCD 3FLOW DATA" and five new subordinate blocks contained within 
the FLOW DATA block which will be headed by: 

BCD 3NETW0RK "Name" 

BCD 3SUBNETW0RK "Name" (Optional) 

BCD 3FLUID LUMP DATA 

BCD 3VALVE DATA (Optional) 

BCD 3FL0W SOURCE DATA 

* 

Table 1 shows the overall organization of the input data blocks including the 
new FLOW DATA block. The five subordinate flow blocks may be input in an arbitrary 
order within the FLOW DATA block. That is, the FLUID LUMP DATA block may be input 
first if desired instead of the order shown. The FLOW DATA block is optional. Thus, 
if the problem being analyzed contains no fluid flow or the fluid flow is being 
handled by another means, there will be no FLOW DATA block. As shown in Table 1, 
the FLOW DATA block will be added between the CONDUCTOR DATA and CONSTANTS DATA 
SINDA input blocks. The flow data input is initiated by the card (starting in 
column 8) 

BCD 3FL0W DATA 
and is terminated by the card 
BCD 3END TLOW DATA 

The NETWORK blocks and SUBNETWORK blocks may be supplied a multiple number 
of times. Each must reference a unique four character name. The user would nor- 
mally supply one NETWORK block for each flow system being analyzed. At least one 
NETWORK block is required if a BCD 3F LOW DATA card exists. The NETWORK blocks may 
or may not reference SUBNETWORK blocks, but if one is referenced the data must be 
supplied in the block referenced by name in the NETWORK data. Division of a network 
into subnetwork elements is sometimes desirable to permit more efficient analysis 
on some problems. The NETWORK DATA, FLUID LUMP DATA, FLOW SOURCE DATA blocks 
are all required when a BCD 3FL0W DATA card exists. Each of these four blocks 
contain the required information for the entire flow oroblem. 

Each of the subordinate flow blocks are discussed in the following sections. 


* 


See Reference 7 for input other than FLOW DATA block 
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TABLE 1 

SINFLO INPUT BLOCKS 


BCD 

END 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 

BCD 

END 

BCD 

END 

I 


BCD 

END 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 

END 

BCD 


3THERMAL LPCS 
3NODE DATA 

3S0URCE DATA (OPTIONAL) 
3CONDUCTOR DATA 


3FLOW DATA * 

3NETWORK (or SUBNETWORK) Name 1 

3NETW0RK (or SUBNETWORK) Name 2* 



— Variable number of NETWORK or 
SUBNETWORK blocks. Must be a 
NETWORK block for each flow 
system. 


3NETW0RK (or SUBNETWORK) Name n* 

3FLUID LUMP DATA 

3VALVE DATA (OPTIONAL) 

3FLOW SOURCE DATA 

3END FLOW DATA 
3CONSTANTS DATA 

3ARRAY DATA 

3EXECUTI0N 

3VARIABLES 1 

3VARIABLES 2 

30UTPUT CALLS 

SEND OF DATA 


> FLOW DATA block optional 

J 


* 

Unique four character name for each NETWORK or SUBNETWORK 
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4.1 


NETWORK and SUBNETWORK Formats 

The fluid flow tubes, pressure nodes connected by the tubes, and 
fluid lump/tube lump pairs contained in each tube are input in the NETWORK 
or SUBNETWORK data blocks. In addition to this connections data, the fluid 
thermo physical property data, network solution parameters, the value of accel- 
eration of gravity and specified pressure nodes and values are input in the 
NETWORK data block. One NETWORK data block must be supplied for each fluid 
system or loop. Connections data for any subnetwork elements are supplied in 
the SUBNETWORK blocks. As many NETWORK and SUBNETWORK data blocks as required 
are supplied. 

The input formats for the NETWORK and SUBNETWORK blocks are shown 
in Table II. Each network or subnetwork block is headed by 
BCD 3NETW0RK Name 
or BCD 3SUBNETW0RK Name 

where the BCD starts in Column 8 and Name is any Alpha/Numeric word up to 
four characters which is different from the name of any other network or 
subnetwork. Each network or subnetwork block is terminated by an END starting 
in Column 8. The data values between the heading and the END card may be input 
free field between columns 12 and 72 consistent with the standard SINDA format. 
Each NETWORK block is a main network for a system. Thus, there are as many NET- 
WORK blocks as there are systems in the problem. The NETWORK block contains (1) 
system data including information referencing fluid property data, system 
solution parameters and specified pressure nodes for the network (the SUB- 
NETWORK blocks do not contain any of this information) and (2) the tube/pressure 
node connections and the fluid and tube temperature lumps in each tube. As 
shown in Table II, the systems data including the property data, solution 
parameters and specified pressure nodes for the network are input immediately 
following the heading card. These items may be input in any order and may be 
supplied one to a line or several to a line. The property values are identified 
by mnanonic names: CP for specific heat, RO for density, MU for viscosity, KT 
for thermal conductivity, and H for enthalpy. The property names are followed by 
an equal sign which is followed by a reference to the property value. For 
example, the property values could be supplied by: 

CP = 0.25, MU = A25 

RO = A37, KT = .073, H = A8 
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TABLE 2 

INPUT FORMAT FOR THE NETWORK AND SUBNETWORK DATA BLOCKS 


BCD 3NETW0RK Name 1 

CP = AXX , RO = AXX, MU = AXX, KT = AXX, GC = XX. XXX, H = AXX 
MPASS = XX, TOL = XX, MXPASS = XX, FRDF = O.XX, KOP = X 
P(N) = XX. X, END 


NTt, NPF], NPT! = FL 11# TL 11# FL 12 , T 12 , 
NT 2 , NPF 2 , NPT 2 = FL 2 i , TL 2 1 , (f 22 , tl 22 , 
NT 3 , NPF 3 , NPT 3 = (FL 31 , TL 3 ! , FL 3 n, TL 3 n 


— F] n , TL] n , END 
FLpns TL2n)» END 
, IFL, 1TL), END 


NT n , NPFp, NPT n “ FL n ], TL n ] — EL nn , TL nn , END 
END 

BCD 3 SUBNETWORK Name 2 

NT], NPF], NPJ 2 » END 

I 1 1 

I t I 

NT n , NPF n , NPT 2 > END 

END 

BCD 3NETW0RK Name 3 

c p - — P ( N ) = XX. X, END 

NT, NPF, NPT = , END 

I 

I 


END 

BCD 3SUBNETWORK Name 4 
END 


The following definitions apply to the above: 


Name i 
CP 
RO 
MU 


any unique four character name 
indicates specific heat value 
indicates density value 
indicates viscosity value 
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TABLE 2 (CONT'D) 


KT 

H 

A XX 
GC 

MPASS 

MXPASS 


TOL 


FRDF 


P(N) 


NT i 

NPFi 

NPTi 

FL ij 

TLij 

IFL 

ITL 

KOP 


indicates thermal conductivity value 
enthalpy value 

array in the ARRAY DATA with actual value of XX 
acceleration of gravity in the desired units. 

Default value = 416962080. 

a pressure/flow solution is performed every MPASS 
temperature iterations. Default value - 1 
maximum number of passes permitted in the balancing 
loop to obtain a pressure/flow solution or any given 
network. Default value = 100 

the solution tolerance on the fraction of change of 
flow rates from one pass in the flow solution to the 
next. Default value = .01 

flow rate damping factor which is a value between 0.5 
and 1.0 to aid the convergence of the flow solution. 
Default value - 0.5 

references the value of the seed fled pressure for 
pressure node N 

tube number i which connects pressure nodes WPF. and 

NPT. 1 

from pressure node number for tube no. i 

to pressure node number for tube no. i 

jth fluid lump in ith tube 

jth tube lump in ith tube 

increment for generating fluid lump numbers 

increment for generating tube lump numbers 

checkout print code (Default value = 0) 

0 : no checkout print is obtained for the network 

1 : a checkout print is obtained for the network 
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The value to the right of the equal may either be constant or reference an 
array in the ARRAY DATA. In the above example., the specific heat and thermal 
conductivity reference constant values of 0.25 and 0.073 while the viscosity 
and density reference arrays 25 and 37 in the array data. The enthalpy curve is 
supplied by array 8. The arrays referenced must be temperature dependent. The 
solution parameters which may be input in the systems data are MPASS, MXPASS, KOP, 

TOL, and FRDF (these are defined in Table II). These items are input by the same 
format as the property data except only integers are permitted for MPASS, KOP and 
MXPASS and only real numbers are permitted for TOL and FRDF. Any or all of the 
five solution parameters may be omitted and default values will be supplied. The 
default values are MPASS = 1, MXPASS * 100, KOP = 0, TOL = .01, FRDF = 0.5, when 
omitted. The acceleration of gravity. GC, is supplied in the systems data. This permits 
the user to analyze the flow problem in any desired units. The default value of GC 
is 416962080, ft/hr . Values of GC for various problem units are given in Table III. 

The specified pressure nodes and their pressure values are also 
supplied in the systems data. For example, if pressure node 34 is set at 
14.7, the input would read 
P(34) = 14.7 

The system data input is terminated by an END similar to ARRAY DATA input. 

An example of the systems data input is 
BCD 3NETW0RK SYSTM1 

CP = 0.25, MU = A25, R0 = 37, KT = 0.073, H = A8 
TOL = 0.01, FRDF = 0.55, MPASS = 2, MXPASS = 120 
GC = 32.173, P ( 34 ) - 14.7, END 
SUBNETWORK input blocks contain no systems data. 

Tube to pressure node connections are supplied in the NETWORK blocks 
following the systems data described above and in the SUBNETWORK blocks. The 
format for the input of the tube cards is: 

NT, NPF, NPT = FL-j , TL-j , FLg, TLg, 

where NT, NPF, NPT are the tube number, "from" pressure node and "to" nodes 
respectively. FL, TL are the fluid lumps/ tube lumps pairs contained 
in the tube. 

A variation of the above format is available for the input of groups of fluid 
and tube lumps with a fixed interval between their numbers. The format for 
this option is: 
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TABLE 3 VALUE OF GC FOR VARIOUS PROBLEM UNITS 
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NT, NPF, NPT = FL-j , TL-, , 


— (FL., 


TL: 


FLj , TL IF, IT),— FL n , TL n , END 


Where the lumps within the parenthesis are being incremented 

FL- is the first fluid lump number of the interval 

FL. is the last fluid lump number of the interval 
J 

IF is the increment between the lump numbers in the interval 
TL^ is the first tube lump number of the interval 
TLj is the last tube lump number of the interval 
IT is the increment between tube lump numbers 
The values of FL.- - FL.. must be a multiple of IF and TL i - TLj must be a mul- 
tiple of IT. If IF and IT are both the integer 1, they may be omitted. 

An example of the input of a tube in the NETWORK block or SUBNETWORK 


is 


8, 3, 5 = (1, 201, 10, 210), END 


This statement indicates that tuoe No. 8 connects pressure nodes 3 and 5 and 
contains temperature fluid lumps 1 thru 10 with adjacent tube temperature lumps 
201 thru 210. A tube card is supplied for each tube in the network. 

Special tube cards are supplied when a subnetwork is referenced from 
a main network. This card consists of a dummy tube number and the first and last 
pressure nodes of the subnetwork on the left of the equal sign and the subnetwork 
name on the right of the equal sign. For example, an input in the main network 
of 

46, 10, 21 = SUB1, END 

references the block with the heading card of 
BCD 3SUBNETW0RK fUBl 

for that portion of the network between pressure nodes 10 and 21. The dummy 
tube number is 46. Subnetwork elements may be referenced from "first level" 
subnetwork elements as well as network elements. 

The input of negative fluid temperature lump numbers on the tube 
cards will indicate that no pressure drop calculations will be made for that fluid 
lump. Negative tube lump numbers indicate no temperature calculations will be made 
on the fluid lumps and tube lumps. This capability is useful for closed loop systems. 
For example, the input for the first tube in a closed system would be 
1, 1, 2 = -200, -297, 97, 297, END 
where fluid lump 200 is the last lump in the system. 
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Each tube must have ai least one fluid lump. This requirement is 
necessary to provide thermal continuity in the network. 

4.2 FLUID LUMP DATA Block Format 

The FLUID LUMP DATA block contains the type data for all fluid 
lumps in all systems. The block is headed by 
BCD 3FLUID LUMP DATA 
and is terminated by 
END 

Where the BCD and END are each in columns 8, 9 and 10 consistent with SINDA 
input convention. The format for this block is the type data for each fluid 
lump type followed by an equal and the fluid lump numbers. The format is 


CSA, WP, FLL, AHT, NHL, MFF, 

FFC, FI, F2 = FLi . FL 2 , — FL n , END 

or 

= (FLi, FL n , INC), END 

or 

= ( FL i , FL n ), END 

or 

= FL], FL 2 , — , (FLi, FLj, INC),— ,FL n , END 


where CSA 
WP 
FLL 
AHT 
NHL 

MFF 


FFC 

FI & F2 
If FI 


cross sectional area to flow 

wetted perimeter 

fluid lump length 

area for convection heat transfer 

a real constant : NHL is the number of head losses 

AXX : XX is an array number of head losses vs Reynolds number 

method for friction factor calculation 

0 : internal calculations used for friction factor 

AXX : XX is an array number of an array of friction factor 

vs Reynolds for Reynolds numbers greater than 2000. 

constant to be multiplied times the friction factor 

FI is a code to determine the method to be used for 

calculating heat transfer 

real number : FI is the convection laminar flow fully 
developed coefficient; F2 is the entry length coefficient 


53 


If FI 


= I : F2 is AXX where XX is an array of Stanton number 

versus Reynolds number from which the convection 
heat transfer coefficient will be obtained 
If FI = 2 : F2 is AXX where XX is an array of heat transfer 

coefficient versus flow rate 
FL -j - ith fluid lump number in the tube 

INC = the increment between lump numbers generated using 

the parenthesis option. If INC is 1 it may be ommited. 

The parenthesis option may be inserted anywhere in the group of fluid nodes 
on the right of the equal sign. That is, lump numbers separated by commas may or 
may not be input before or after the lump generated by the paraenthesis option. 

Also, any number of parenthesis options may be used on one type card. 

The values for MFF, FFC, FI, and F2 may be left off the type cards 
if the default values are desired for al 1 these items. The default values 
are MFF = 0 and FFC, FI and F2 = 1.0. The input would then be: 

CSA, WP, FLL, AHT, NHL = FL -j , FL 2 , — (FLj, FL j , INC), — FL n , END 

4.3 VALVE DATA Input Block (Optional) 

The VALVE DATA input block contains the valve data for all valves 
in all flow systems. (That is, there is only one VALVE DATA block in the 
problem.) The block is headed by the card 
BCD 3VALVF. DATA 
and is terminated by 
END 

Where the BCD and the END cards are each in columns 8, 9, and 10. 

Three types of valves are available to the user: .rate limited, polynomial, 
and switching {see Section 2.2.3). The input required for these valves is: 

Rate Limited 

NV, NTS1 , NTS2 = XI, MODE, XMIN1, XMAXl, E, TSEN1 , TSEN2, DB, RF, RL, END 
Polynomia l 

NV, NTS1 , NTS2 = XI, MODE, XMIN1 , XMAXl, E, TSEN1 , TSEN2, CO, Cl, C2, C3 , C5, VTC, END 
Switching 

NV, NTS 1 , NTS2 = XI, MODE, XMINl, XMAXl, E, NSEN, T1 , T2, END 


5 * 


NV 

NTS1 

NTS2 

XI 

MODE 

XMIN1 

XMAXI 

E 


TSEN1 


TSEN2 


CO, Cl 


DB 

RF 

RL 




The following definitions apply in the above formats: 

Valve number 

Tube number connected to side 1 of the valve 
Tube number connected to side 2 of the valve 
Initial valve position 

Operating mode : 1 - operating; 0 - not operating 
Side 1 minimum position; side 2 maximum position is 
(1.0 - XMIN1) 

Side 1 maximum position; side 2 minimum position is 
(1.0 - XMAXI) 

The valve geometric factor relating pressure drop through 
the valve by 

p 

AP = E (flowrate/valve position) 

Sensor lump for side 1 or set point for side 2; if TSEN1 
is an integer, it identifies the side 1 sensor lump to 
be controlled to (a) the set point for side 1 or (b) the 
sensor lump for side 2 (TSEN2). If the variable is input 
as a floating point number it represents a set point to 
which the side 2 sensor lump will be controlled. 

Sensor lump for side 2 or set point for side 1; if TSEN2 
is an integer, it identifies the side 2 sensor lump to be 
controlled to (a) the set point for side 2 or (b) the sensor 
lump for side 1 (TSEN1). If the variable is input as a float- 
ing point number it represents a set point to which the side 
1 sensor lump will be controlled. 

C2,C3,C4,C5 - Polynomial curve fit coefficients for a curve fit of 

the steady state valve position vs sensed temperature error 
for side 1 : 

X1SS = CO + Cl -AT + C2-AT 2 + C3 • AT 3 + C4-AT 4 + C5-AT 5 
Dead band for the rate limited valve, degrees of temperature 
(Figure 4) 

Rate factor, the rate of change of valve velocity to sensed 

temperature error (dx/d(AT)) as shown on Figure 4 

Rate limit, the maximum valve velocity, Xmax (see Figure 4) 


VTC - Valve time constant as described in Section 2. 2.3. 2. 


If a valve is desired with no time lag, a time constant 
which is very small compared to the problem time increment 
should be input. (VTC must be greater than zero). 


NSEN 

- 

Sensor lump for switching 

valve 


T1 

- 

Side 1 off temperature or 
valve 

side 2 

on temperature for switching 

T2 

“ 

Side 2 off temperature or 
valve 

side 1 

on temperature for switching 

4.4 

FLOW SOURCE Data Block 




The FLOW SOURCE data block contains specification of flow rate for 
all the systems in the problem. The heading card for this block is 
BCD 3FL0W SOURCE DATA 
and is terminated by 
END 

Three types of flow specifications are available. These are: 

(1) flow rate as a function of time; (2) pressure rise as a function of flow 
rate specified by a tabulated curve; and (3) pressure rise as a polynomial 
function of flow rate. The input for each of these is given below. 

Flow as A Function of Time 
NPI , AW, END 

Pressure Rise as A Tabulated Function of System Flow Rate 
NPI, NPO, ADP, END 

Pressure Rise as A Polynomial Function of Flow Rate 
NPI, NPO, CO, Cl, C2, C3, C4, END 

where: 


NPI 

AW 


NPO 

ADP 


system inlet pressure node 

W is an array number of an array which gives tabulated 

flow rate vs time if input as AXX 

AW is constant imposed flowrate for node NPI if AW 

is a floating point number. 

system outlet pressure node 

DP is an array number of an array which gives tabulated 
pump pressure rise as a function of flow rate 


CO, Cl ,C2,C3,C4 


= polynomial curve fit constants for pressure rise 
as a function of flowrate, i.e., 

AP ~ CO + Cl *w + C2'W^ + C3*fii^ + C4*ti/^ 

The value of AW may be input as a floating point number if a constant system 
flow rate is desired. 

Example of Flow Input 

An example of the flow input described in Sections 4.1 and 4.4 is 
in Table IV. This table gives the flow input for the sample problem 
in Section 6,0 . 


4.5 


given 

given 


I 

I 

I 


1 

I 

B 

I 

I 
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TABLE 4 

FLOW DATA INPUT FOR SAMPLE PROBLEM 


rcd in nw data 

BCn 3NETW0KK MAIN 

GC = q . 1 73 1 2Eb , C P 55 A 1 « RO **2. MU»A5. K T * A 6 » MPA55»1» H*A8 

T0L««0l»MXPASSsl3QtFRDF’*»7,Pt24!*0.»END 


37 * 

23 * 

Z8 

SB 

1 l / .3 17 

6 

FND 

. 

2 . 

2 Q 

3 

= 

98,29a 

ft 

FND 


3 p 

3 , 

9 

= 

99 .299 

• 

fnd 

i 

j 

3B» 

9 » 

1 7 

a 

SUB l 

• 

END 


7 6 , 

i 7 , 

l 8 

3 

1 Q6 .306 

ft 

fnd 

1 1 i 

3. 

9 

3 

107 , 307 

0 

fnd 

i 

i 

39, 

9 » 

22 

B 

SUB2 

ft 

fnd 

i 

39 t 

22 • 

18 

* 

llHOH 

• 

FND 

35, 

Id, 

23 

a 

1 15,315 

ft 

END 

; 

3 6 • 

2 * 

23 

s 

116,316 

ft 

END 

1 

1 * 

t i 

2 

3 

-200.-297 ,97,297 

ft 

FND 

i 


FND 

PCD 35 UR NETWORK SURl 


*1 . 

8 * 

6 

s 


1 00 , 30n 

0 

END 

5 o 

6 * 

7 

3 

i 

7 , 207 , 

12.212), 

FND 

6 , 

6 » 

7 

3 

< 

1 » 2^ 1 , 

6,206) , 

fnd 

7. 

8 6 

5 

3 

i 

13,213, 

lb, 218) , 

FND 

8. 

8 , 

5 

3 

( 

19,219, 

28 , 228 ) , 

FNI) 

9 . 

5 • 

7 

3 


101.301 

» 

fnd 

10, 

7, 

8 

JX 


102,302 

« 

FND 

19, 

8, 

1 8 

X 


103,303 

* 

FND 

20 , 

18, 

16 

X 


1 08 , 308 

» 

FND 

21 , 

16, 

l 7 

a 

( 

37 ,237 , 

92,29?), 

END 

22, 

16, 

1 7 

a 

t 

83,283, 

9b • 29ft ) • 

END 

23. 

18, 

15 

B 

( 

31,231, 

36,236). 

fnd 

28, 

18, 

15 

a 

i 

25,225, 

30,230) » 

END 

25, 

15, 

1 7 

s 


I 05,305 

» 

END 


FNO 

RCO 3 SIJRNETWORK SURZ 


1 2 . 

9 , 

1 1 

S 


t Ob » 308 

ft 

END 

13 , 

1 1 . 

1 2 

c 

1 

6 1 , 2 6 | , 

66 ,266 > . 

END 

1 9 , 

1 1 , 

1 2 

E 

t 

6 7 ,2 6 7 , 

7 2 » 2 7 2 ) * 

END 

15 , 

9 , 

10 

X 

( 

55,255 , 

60 , 260 ) , 

end 

1 6 » 

9 . 

10 

X 

( 

99 , 299 , 

59 , 259 ) . 

FND 

1 7 , 

10 , 

1 1 
* M- 

c 


1 09 ,309 

ft 

END 

1 8 , 

1 2 , 

1 3 

X 


110,310 

ft 

END 

27 , 

13 . 

19 

n 


111,311 

ft 

END 

28 , 

19 , 

2 l 

X 


1 U .312 

ft 

END 

29 . 

2 1 , 

22 

X 

l 

79 , 279 , 

84 *284 ) * 

END 

30 , 

2 1 , 

22 

X 

I 

73 , 2 / 3 , 

78 , 278 ) , 

FND 

3 l . 

1 9 „ 

20 

X 

( 

85 , 285 . 

90 * 290 ) , 

FND 

32 . 

19 , 

20 

X 

( 

91 , 291 , 

96 , 294 ) , 

END 

33 , 

20 , 

27 

G 


113,313 

• 

FND 


END 


} I 
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TABLE 4 (CONT'D) 


BCD 3FLUID LUMP DATA 


o.ooioob, 

0, 

1125, 12.0 , 

1«35 « 

0,0 ® 







1 

t 6 

* 19 

. 29 * 

25 t 

3n 


93 

, 

98 


<49 

. 59 

» 67 

, 72 * 

73 » 

78 


91 

, 

96 

n.oon93fi , 

n , 

36 , 

3*25, 

1.17 , 

1 17.0 

8 





( 

2 

, 51 

t 1 8 

* lilt 

19 , 

17) 

, t 

20 

, 

23) 

{ 

76 

« 29 ) 

, f 32 

i 35 ) » 

38 • 

91 > 
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USER SUBROUTINES 


This section describes the user subroutines which have been developed 

and modified by v$D for SINDA. Table 5 summarizes the subroutines and the page 

that each description of usage is found. 

The subroutine inputs rely upon the ability to convert from actual 

array, node, and conductor numbers to relative numbers in the array data. To 

use the capability the user may supply an actual array number, node number, or 
conductor number by preceding the actual number with *A, *T, or *G respectively . 

This causes the preprocessor to repine the entry with the relative number. Con- 
sider the example for array number 2 shown below 

2, *A14, *T5, *G7, END 

In this example, following the preprocessor phase, *A14 would be replaced by the 
location in the A array of the array number 14 data, *T5 would be replaced by 
the relative node number for actual node number 5, and *G7 would be replaced by 
the relative conductor number for actual conductor number 7. 

In addition, revisions have been made to some of the temperature solution 
subroutines to interface with the Fluid Hybrid solution subroutines. The following 
subroutines were revised: 


CINDSL 

Steady State 

CINDSS 

Steady State 

CNBACK 

Backward Differencing 

CNFWBK 

Mid-differencing 

CNFRWD 

Forward Differencing 

CNFAST 

Forward Differencing 

FWDBCK 

Mid-differencing 

SNDSNR 

Steady State 

SNFRDL 

Forward Differencing 

SNFRWD 

Forward Differencing 

STDSTL 

Steady State 
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TIMCHK 1 01 

WPRINT . 102 


61 


SUBROUTINE NAME: 


ACOMB 


See description for usage of OCOMB. 



SUBROUTINE NAME: 


CABIN 


PURPOSE : 

This subroutine performs a thermal and mass balance on a cabin air system. 
The cabin air is assumed to be a two component gas mixture with one con- 
densible component and one noncondensible component. The cabin air is 
assumed to be well mixed so that the temperature and specific humidity are 
constant throughout. The cabin may contain any number of entering streams 
each with different temperature and humidity conditions. The cabin air 
may transfer heat to any number of nodes in its surroundings with the heat 
transfer coefficient obtained by one of the the three options: 

1. User input coefficient 

2. Relations for flow over a flat plot 

3. Relations for flow over a tube bundle 

The relations describing the second and third options are given in Section 
2.1.7. The mass transfer coefficient for determining the rate of condensa- 
tion or evaporation is determined by the Lewis relation which relates the 
mass transfer coefficient directly to the convection heat transfer coeffi- 
cient. By the Lewis Relation, if the diffusion coefficient is approximately 
equal to the thermal diffusi vity , the Sherwood number is approximately equal 
to the Nusselt number, thus giving a direct relation. (See Section 2.1.7 
for details). Mass and heat transfer rates are determined at each node 
that interfaces the cabin gas as well at entering and exiting streams and 
a new cabin gas temperature and humidity is determined each iteration based 
upon the heat and mass balance. An account is kept of the condensate on the 
walls when condensation occurs but the condr-"-ate is assumed to remain 
stationary and not flow to other wall no 

Limits are applied when necessary to prevent more condensation than the 
vapor existing under severe transient condition and to prevent evaporation 
of more liquid than exists at each wall lump. 

As many cabins as desired may be analyzed in a given problem, but each must 
contain separate input information. 

RESTRICTIONS : • 

CABIN must be called in VARIABLES 1. 

CALLING SEQUENCE : 

CABIN(A(IC) TC, TC, K1 , K2) 

The following definitions apply to the above calling sequence: 

A is an array containing arrays numbers which 

contain cabin input information 
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TC 


The cabin gas temperature which must be a 
boundary node 

Storage locations needed by CABIN 


K1 >K2 

The A array has the following format where the *A procedure is used: 
A(IC) , IF s PR,CNjH s FP,TB,SP,END 


Where IF 


PR 


CN 


H 


Identifies an array containing the entering flow 
rate information. The format of the array is: 

IF(IC) # NS,FR 1 ,PSI,TE 1 ,FR 2 ,PSI 2 JE 2 - — FR ns >PSI ns ,TE ns 

Identifies an array identifying array numbers 
for property values. The format of the array 
is: 

PR( I C) ,NFLC jNMUO, NMUV ,NCP0, NCPV ,NK0 ,NKV ,NLAT 

Identifes an array containing pertinent constants. 
The format of the array is: 

CNttCl jRA,RV,VC,PC ,XC,WV .PSIC,P0 jT0,C0NV 

Identifies an array containing node numbers 
and convection heat transfer coefficient 
values for nodes surrounding the cabin gas. 

The format of the array is: 

H(LC) s LNi > HAi s LN 2 , HA 2 , LN n i , HA nl 


FP Identifies an array containing node numbers and 

information to permit calculation of convection 
coefficients for flat plates. The format is: 

FP(T.C) ,LN ] ,XX 1 .XI, »AI, .VIW0, 5 LN 2 ,XX 2 ,XI 2S AI 2 , 

VIN 0 2 , LN n2 .XX n2> XI n2 .AI n2 ,VIW0 n2 

TB ' 'Identifies an array containing node numbers and 

information to permit calculation of convection 
coefficients for tube bundles. The format is: 

TB(IC) .LM, ,DI AI, .VIW0-J ,LN 2 ,DI 2 ,AI 2 ,VIW0 25 LN n3> 

DI n3 ,AI n3’ VIWO n3 

SP Identifies an array which contains working space 

equal to or greater than three times the sum 
of the number of nodes with input heat transfer 
coefficients plus the number using flat plot 
relations plus the number using tube bundles. 

The following symbol definitions apply in the above: 


F4 


NS 


FR i 


PSI. 


TE. 


NFLC 

NMUO 

NMUV 

NCPO 

NCPV 

NKO 

NKV 

NLAT 

RA 

RV 

VC 

PC 

XC 

WV 

PSIC 



HA 


nl 

n2 


Number of incoming streams 

Entering flow rate for stream i 

Specific humidity for entering stream i 

Temperature of entering stream i 

Curve number for circulation flow rate vs time 

Curve number for noncondensible viscosity vs 
temperature 

Curve number for condensible viscosity vs 
temperature 

Curve, number for noncondensible specific heat 
vs temperature 

Curve number for condensible specific heat vs 
temperature 

Curve number for noncondensible thermal conduc- 
tion vs temperature 

Curve number for condensible thermal conduction 
vs temperature 

Curve number for latent heat of condensible vs 
temperature 

Gas constant for non-condensible component 
Gas constant for condensible component 
Cabin volume 
Cabin Pressure 

Molecular weight ratio, Mv/Mo 
Initial vapor weight in cabin 
Initial specific humidity for cabin 
Cabin wall lump 

Heat transfer coefficient times area 

Number of wall lumps which have input HA values 

Number of wall lumps which have HA calculated 
by flat plate relations 

Number of wall lumps which have HA calculated 
by tube bundle relations 
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n3 


Distance from leading edge for flat plate 
heating for ith flat plate node 

Length of flat plate in flow direction for ith 
flat plate node 

Heat transfer area for flat plate or tube node 

Tube outside diameter for tubes in the bundle 
for ith tube node 

Ratio of velocity at the lump to the circulation 
flow rate 

The reference temperature to be used for esti- 
mating the saturation pressure of the condensi- 
ble component. Should be near the range of 
saturation temperature expected 

The saturation pressure at To for the condensible 
component 

Conversion factor to make the quantity XLAM/Rv/To 
dimensionless where XLAM is the latent heat of 
vaporization and Rv is the gas constant for the 
vapor. If XLAM is BTU/1 b » Rv is FT-LB/°R and 
To is °R, C0NV=778. 


SUBROUTINE NAME: 


CRVIMT 


PURPOSE : 

This subroutine performs an integration of the doublet array, A, and 
stores the results in doublet array B. The independent variables (the odd data 
valves) of the A array are transferred directly to the B array. The dependent 
variables of the B array are calculated by 

B ( 2) = 0.0 

B(2*N) = B(2*(N-1 ) ) + G.5*[A(2*N) + A(2*(N-1)] 
*[A(2*N-1)-A(2*(N-1)-1] 

N =2, NP 

where NP = number of points in the A array 
(half the integer count) 

This subroutine was written primarily for integration of specific heat arrays 
to obtain enthalpy arrays but could be used for integration of any dependent 
variable over the independent variable range. 

RESTRICTIONS : 

Space in B array must be exactly equal to the space in the A array 
Must be at least two points in A array (i.e., the integer count must be greawr 
than or equal to 4). 

CALLING SEQUENCE : CRVINT(A(IC) , B( IC) ) 
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SUBROUTINE NAME: 


CYCLE 


PURPOSE : 

Subroutine CYCLE will automatically extend the range of independent 
variables in either direction for cyclic curves by adding (or subtracting) the cycle 
period to each independent variable when the curve range is exceeded. The total 
input range of the independent variables is assumed to be one period. CYCLE should 
be called prior to interpolation so that the necessary changes may be performed to 
the independent variables. 

RESTRICTIONS : 

None 

CALLING SEQUENCE : CYCLE(X,A,NAME) 

where X - value of independent variable to be used in interpolation 

of A doublet array 

A - doublet array assumed to be cyclic 
NAME - one word Hollerith identifier 
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SUBROUTINE NAME: 


FLOSOL 


PURPOSE : 

Subroutine FLOSOL determines the flow distribution in a set 
of general para 1 lei /series fluid flow tubes so that the pressure drop values 
between any parallel flow paths are equal and flow is conserved. The following 
effects are included in the pressure drop calculations: 

(1) pipe flow friction 

(2) orifices and fittings 

(3) valves 

The effect of temperature dependent properties are included in the calculations. 

The properties are evaluated at the temperature of each fluid lump in each tube 
in evaluating the flow resistance when setting up the equations to be solved. 

A balance is made between the flow/pressure drop characteristics of the system 
and the flow/pressure rise of a pump for each system concurrent with the system 
pressure flow solution to obtain the incoming system flowrate. A detailed dis- 
cussion of the equations and techniques used are described in Section 2.1. General 
flow charts of FLOSOL and supporting subroutines are shown in Figures 6,7, and 8. 

RESTRICTIONS : 

FLOSOL should be called from EXECUTION prior to temperature 
solution call and from VARIABLES 2 for transient problems. For steady state solu- 
tions FLOSOL should be called from VARIABLES 1 and DTIMEU should be set in the 
CONSTANTS DATA if valve operation is required. The system of units used for the 
thermal and flow problem should be consistent. 

CALLING SEQUENCE : FLOSOL 

DYNAMIC STORAGE REQUIREMEN TS: 

Dynamic storage required for FLOSOL is 1 / 2 (NPRN + 

7*NPRN + 12), where NPRN is the maximum of the number of pressure nodes in any 
network. 
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CALCULATE 

VALVE 

RESISTANCES 


SET UP SPACE 
FOR MATRIX, 
RIGHT HAND 
SIDE, ETC. 


ASSEMBLE 

COEFFICIENT 

MATRIX 


STORE 

IMPOSED FLOW 
RATES IN RHS 


STORE 
TNLET FLOW 
RATE IN RHS 


REDUCE 
MATRIX AND 
ADJUST RHS 


CALL 

SYMSOL 

FOR 

PRESSURES 


FLOW 


CALCULATE 

NEW 

LOW RATES 


RETURN 


FI SURE ** FLOW CHAR”* OF FLBAL 






72 








SUBROUTINE NAME : 
PURPOSE: 


FLOTMP 


Subroutine FLOTMP will read the node temperatures, flowrates, 
pressures and valve positions at time TMPTIM from the history tape assigned to Unit U 
generated by subroutine HSTFLO for a previous run on Unit T to initiate a problem at 
tnese conditions. The time to read the tape, TMPTIM is the argument. The subroutine 
should be called in the execution block prior to the call to the temperature solution 
subroutine. 

RESTRICTIONS : 

Must be called in the EXECUTION block prior to the call to the 
appropriate temperature solution subroutine. The history tape must be assigned on 
Unit U. 

CALLING SEQUENCE : 

FLOTMP (TMPTIM) 


73 


SUBROUTINE NAME: 


FLPRNT 


PURPOSE : 

Subroutine FLPRNT will write the values of the DATA array of 
real numbers at 10 to a line. The array is labeled by the variable input HEAD 
which contains 9 six character alpha numeric words. The array location of every 
tenth value in the array is identified to the right of the appropriate line. 


RESTRICTIONS : 

Should be called from OUTPUT. The array must be real. 

CALLING SEQUENCE : 

FLPRNT(DATA(IC) , HEAD(DV) ) 
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SUBROUTINE NAME : FLUX 

PURPOSE : 

Subroutine FLUX permits doublet time variant curve values stored 
on magnetic tape unit NFLXTP to be read into NCRV arrays starting at array 
DATA when the mission time exceeds DQTIME. The flux tape must be generated 
prior to the run using a GE routine LTVFTP. This routine generates the flux 
tape in the following format: 

Record No. 1 
First Read Time 

Record No. 2 

Number of points on first curve (Integer), first curve indeDendent variables, 
first curve dependent variables, number of points on second curve, second 
curve independent variables, second curve dependent variables, etc. for all 
curves. 

Record No. 3 
Second Read Time 

Record No. 4 

Same as Record No. 2 except with new values 

Record No. 5 
Third Read Time 

Etc. until all blocks of data are on tape. 

Subroutine FLUX writes the values from the appropriate NFLXTP record into 
the arrays defined by DATA and NCRV in the proper doublet array format. 

Flux values should be input into the heat flux arrays (DATA] — DATA^Q^y) 
initially if the user doesn't want the values to be read from the tape 
at the start of the problem. The value of QTIME should initially be the 
value of the time the first read is desired. 

RESTRICTIONS : 

The following restrictions apply: 

(1) The initial block of curve data must be input on cards or data 
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CALLING SEQUENCE : 

FLUX ( NFLXTP , DATA , NCRV. IXJTlMt , QUMO 

v/here 


NFLXTP 

DATA 

NCRV 

DQTIME 

QTIME 


logical unit to which the r 1 u v U|w> i\ assiuntid Mtn> I he 

supplied by a user constant. 

starting location (IC) for Mux curves 

number of flux curves to he updated hum Hie Mu* lope 

time scale shift for flux curves l)(|IIMl Is added In timli 

independent value for each I lux cur ve r ead hum IIM Ml* 

the last point on the latest set of flux turves I 

NFLXTP. (QTIME = FLXT1M t IKJIIMI , where IIM1H is n,c 

read from the flux tape) must he supplied \>/ use* -mr. i<mi 


SUBROUTINE NAMES: 


GENOUT, GENI OR GENR 


PURPOSE : 

These subroutines print out arrays of numbers 10 to a line. 
GENOUT prints either real numbers, integer or both. GENI and GENR print 
integers and real number arrays respectively. The integers are written in an 
19 format and the real numbers in an E12.4 format. 

RESTRICTIONS : 

GENI writes arrays of integers only. GENR writes arrays of 


real numbers only. 
CALLING SEQUENCE: 



GENOUT. (A, ISTRT, ISTP, 'NAME*) 


GENI 

(A, ISTRT, ISTP, 'NAME*) 


GENR 

(A, ISTRT, ISTP, 'NAME') 

where A 

- 

is the array location 

ISTRT 

- 

is the first value in A being written 

ISTP 

- 

is the last value in A being written 

‘NAME' 

- 

is a title of 22 Hollerith words for identification 
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SUBROUTINE NAME : HSTFLO 

PURPOSE : 

Subroutine HSTFLO stores the problem time, the pressures of 
all pressure nodes, the valve positions for all valves, the flowrates for all tubes, 
and the temperatures of all temperature nodes at an input interval on a magnetic tape 
(the history tape) mounted on Unit T. The number of records written on the history 
tape is the number of history intervals plus two. The first record contains a title, 
an integer count of the number of items to be written for each of the five categories 
(pressure drops, pressures, valve positions, flowrates, and temperatures), the actual 
tube numbers, actual pressure node numbers, actual valve numbers, and the actual 
node numbers in order of relative numbers. The second through the next-to-last records 
contain the history records with one for each time point and the last record is the 
same as the next-to-last except the time is negative. The argument to HSTFLO is the 
history tape writing interval, TINC. 

The format for the history tape is as follows: 

Record No. 1 

Tital (written internally) is 12A6 format 0,0, 0,0,0, number of tubes, number of pressure 
nodes, number of valve positions, 0,0,0, number of tubes, 0,0, number of nodes, actual 
tube numbers in increasing order, actual pressure node numbers in increasing order, 
actual valve numbers in increasing order, and actual node numbers in increasing order 
of relative node numbers. 

Record No. 2 

Initial problem time, pressure drops, pressures, valve positions, flowrates, node 
temperatures 

Record No. 3 

Second history time, pressure drops, pressures, valve positions, flowrates, node 
temperatures 


I 


Record No. N+l (Where N = number of history time slices to be written) 

Last history time, pressure drops, pressures, valve positions, flowrates, node temperatures 
Record No. N+2 

Same as last record except time is negative 
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RESTRICTIONS : 

Should be called in VARIABLES 2. An output history tape should 
be mounted on unit T. Subroutine TMCHK must be in VARIABLES 2 prior to the call to 
Subroutine HSTFLO if TIMCHK is ca^ed in the problem. 

If the backup feature is used in VARIABLES 2, the call to subroutine 
HSTFLO should not be made until the last pass to avoid nonincreasing time records 
or invalid data. For example: 

BCD 3VARIABLES 2 


F IF {T (16) .LT. TMAX) BACKUP = 1. 


F IF {BACKUP .GT. C, GO TO 10 

HSTFLO (.01) 

F 10 CONTINUE 
END 

CALLING SEQUENCE : 

HSTFLO(TINC) 


79 


SUBROUTINE NAME: 


HXEFF 


PURPOSE : 

This subroutine obtains the neat exchanger effectiveness either from a 
user constant or from a bivariant curve of effectiveness versus the flow rates on 
the two sides. The effectiveness thus obtained is used with the supplied flow rates, 
inlet temperatures and fluid properties to calculate the outlet temperatures using 
the methods described in Section 2. 1.6. 4. The user may specify a constant effec- 
tiveness by supplying a real number or may reference an array number to specify the 
effectiveness as a bivariant function of the two flow rates. The user supplies flow 
rates, specific heat values, inlet temperatures and a location for the outlet 
temperatures for each of the two sides. The flow rate an-ay may be referenced to 
obtain flow rates and the temperature array may be used for temperatures. The 
specific heat values may be supplied as a temperature dependent curve or a constant 
value may be supplied. The user also identifies enthalpy curves for each side 
which may be generated from the specific heat curve with user subroutine CRVINT. 


RESTRICTIONS : 

HXEFF should be called in the VARIABLES 1 block. The value for EFF, 
the first argument must never be zero. T 0(Jt ^ and T ^ must be boundary nodes. 


CALLING SEQUENl * HXEFF ( EFF, W1 ,W2,CP1 ,CP2,TIN1 ,TIN2,T0UT1 ,T0UT2,H1 ,H2) 

where EFF - is (1) the effectiveness if real, (2) a curve number of 

a bivariant curve of effectiveness versus W1 and W2 if an 
array 

W1 ,W2 - are the flow rates fot side 1 and 2 respectively. May 

reference the flow rate array, W {I)where I is the tube 
number 
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CPI 9 CP2 - 

TIN1 ,TIN2 - 
TOUT!, TOUT 2 - 

H1.H2 






^Ym"r rr 




are the specific heat value for side 1 and side 2 fluid 
respectively. Constant values may be input or arrays 
may be used for temperature dependent properties, 
are inlet lump temperatures - Usually T(IN1) and T(IN2) 
where INI and IN2 are the inlet lumps on side 1 and side 2 
are the outlet lump temperature locations sides 1 and 2 
where the calculated values will be stored. Must be 
boundary nodes. 

are arrays which give enthalpy vs temperature for sides 1 
and 2 respectively. 
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SUBROUTINE NAME: 


HXCNT 


PURPOSE : 

This subroutine calculates the heat exchanger effectiveness using the 
relation described in Section 2. 1.6.1, for a counter flow type exchanger. 

The value of UA used in the calculations may be specified as a constant 
by supplying a real number or it may be specified as a bivariant function of the two 
flow rates by referencing an array number. The user supplies flow rates, specific 
heat values, inlet temperatures and a location for the outlet temperatures for each 
of the two sides. The flow rate array may be referenced to obtain flow rates and 
the temperature array may be used for temperatures. The specific heat values may be 


supplied as a temperature dependent curve or a constant value may be supplied. 
The user also identifies enthalpy curves for each side which may be generated 
from the specific heat curve with user subroutine CRVINT, 


RESTRICTIONS : 

HXCNT should be called in the VARIABLES 1 block. The value of UA, the 
first argument, must never be zero. T QUt ^ and T QUt2 must be boundary nodes. 

CAL l I NS SEQUENCE : HXCNT(UA,W1 ,W2,CPl ,CP2,TIN1 ,TIN2,T0UT1 ,T0UT2,H1 ,H2) 

where UA is (1) the heat exchanger conductance if real, (2) a 

curve number of a bivariant curve of conductance versus 
W1 and W2 if an array 

W1,W2 are the flowrates for side 1 and side 2 respectively. 

May reference the flowrate array ,W (I) where is the 

tube number. 

* 

CP1,CP2 are the specific heat values for side 1 and 2 fluid 

respectively. Constant values may be input or arrays 
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TOUT! 


HI ,H2 


may be used for temperature dependent properties. 

TGUT2 are the outlet lump temperature locations (sides 1 and 
2) where the calculated values will be stored. Must be 
boundary nodes. 

are arrays which give enthalpy vs temperature for sides 
1 and 2 respectively. 
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SUBROUTINE NAME : HXCOND 

PURPOSE : 

This subroutine performs thermal analysis on a condensing 
heat exchanger using relations described in section 2. 1.6. 5. The effectiveness 
may either be supplied as a constant or as a trivariant function of humidity* 
flow rate of the gas, and flow rate of the coolant. CRVINT may be used to 
integrate the specific heat curves to produce the enthalpy curves. 

RESTRICTIONS : 

HXCOND should be called in the VARIABLES 1 block. The value 
for EFF, the first argument, must never be zero. TGOUT, and TCONI must be 
boundary nodes. 

CALLING SEQUENCE : HXCOND ( EFF ,WG ,WC ,NHG , NHC ,TG IN ,TC I N , PS II N , P , XLAM ,XMIMO , 

PSI0UT,WL ,TG0UT,TC0UT) 

v.nere EFF is (1) the effectiveness if real, (2) a curve number of a 

tri variant curve of effectiveness versus PSIIN, WG, and WC 
WG is the flow rate of the gas 

WC is the flow rate of the coolant 

NHG is the enthalpy curve for the gas 

NHC is the enthalpy curve for the coolant 

TGIN is the temperature of the incoming gas 

TCIN is the temperature of the incoming coolant 

PSIIN is the humidity of the incoming gas 

P is the total gas pressure 

XLAM is the latent heat of vaporization 

XMDMD is the molecular weight ration M y /M Q 

PSlOUT is the outlet humidity 

t 
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WL 

TGOLIT 

TCOUT 


is the flow rate of the liquid 

is the temperature of the outgoing gas 

is the temperature of the outgoing coolant 
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SUBROUTINE NAME: 


HXCROS 


PURPOSE : 

This subroutine calculates the heat exchanger effectiveness using the 
relations described in Section 2. 1.6. 3, for a cross flow type ex- 
changer. The value of UA used in the calculations may be specified as a constant 
by supplying a real number or it. may be specified as a bivariant function of the 

two flow rates by referencing an array number. Any one of the following four 
types of cross flow exchangers may be analyzed. 

1) Both streams unmixed 

2) Both streams mixed 

3) Stream with smallest MCp product unmixed 

4) Stream with largest MCp product unmixed 

The type is specified by the last argument in the call statement. The user 
supplies flow rates, specific heat values, inlet temperatures and a location for 
the outlet temperatures for both sides. The flow rate array may be referenced 
to obtain flow rates and the temperature array may be used for temperatures. The 
specific heat values may be supplied as a temperature dependent curve or a constant 
value may be supplied. The user also identifies enthalpy curves for each side 
which may be generated from the specific heat curve with user subroutine CRVINT. 

RESTRICTIONS : 

HXCROS should be called in the VARIABLES 1 block. The value for UA, the 
first argument, must neve" be zero. T Q(J ^ and T Qut2 must be boundary nodes. 

CALLING SEQUENCE : HXCR0S(UA,W1 ,W2,CP1 ,CP2,TIN1 ,TIN2,T0UT1 ,T0UT2,K,H1 ,H2) 

where UA is (1) the heat exchanger conductance if real, (2) a 

curve number of a bivariant curve of conductance versus 
W1 and W2 if an array. 
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W1 ,W2 are the flow rates for side 1 and 2 respectively. Pay 

reference the flow rate array, W (I) where I is the 
tube number. 

CPI ,CP2 are the specific heat values for side 1 and side 2 fluid 

respectively. Constant values may be input or arrays may 
be used for temperature dependent properties 

TIN1,TIN2 are inlet lump temperatures - Usually T(IN1) and 

T(IN2) where INI and IN2 are the inlet lumps on 
side 1 and side 2 

TOUT1 S TOUT2 are the outlet lump temperature locations (sides 1 and 
2) where the calculated values will be stored. Must 
be boundary nodes 

K is the code specifying type of cross flow exchanger: 

Both streams unmixed : K = 1 
Bo ti» streams mixed : K = 2 
Stream with small WCp unmixed : K = 3 
Stream with large WCp urmixed : K = 4 

HI ,H2 are arrays which give enthalpy vs temperature for sides 

1 and 2 respectively 
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SUBROUTINE NAME : HXPAR 

PURPOSE : 

This subroutine calculates the heat exchanger effectiveness using the 
relations described in Section 2. 1.6. 2, for a parallel flow type exchanger. 

The value of UA used in the calculations may be specified as a constant 
by supplying a real number or it may be specified as a bivariant function of the two 
flow rates by referencing an array. The user supplies flow rates, specific heat 
values, inlet temperatures and a location for the outlet temperatures for each of 
the two sides. The flow rate array may be referenced to obtain flow rates and the 
temperature array may be used for temperatures. The specific heat values may be 
supplied as a temperature dependent curve or a constant value may be supplied. 

The user also identifies enthalpy curves for each side which may be generated 
from the specific heat curve with user subroutine CRVINT. 


RESTRICTIONS : 

HXPAR should be called in the VARIABLES 1 block. The value for UA, the 
first argument, must never be zero. T ^and T Qut2 must be boundary temperatures. 

CALLING SEQUENCE : HXPAR{UA,W1 ,W2,CP1 ,CP2,TIN1 ,TIN2,T0UT1 ,T0UT2,H, ,H2) 

where UA is (1) the heat exchanger conductance if real, (2) a curve 

number of a bivariant ^urve of conductance versus W1 
and W2 if an array. 

W1,W2 are the flow rates for side 1 and 2 respectively. May 

reference the flow rate array, W {l)where I is the tube 
number 

CPI ,CP2 are the specific heat values for side 1 and side 2 fluid 

respectively. Constant values may be input or arrays may 
be used for temperatures dependent curves. 


TIN1.TIN2 


are inlet lump temperatures - Usually T(IN1) and T(IN2) 
where INI and IN2 are the inlet lumps on side 1 and side 2 
T0UT1,T0UT2 are the outlet lump temperature locations (sides 1 and 2) 
where the calculated values will be stored (should be 
boundary temperatures) 


HI ,H2 


are arrays which give enthalpy vs temperature for sides 1 
and 2 respectively 


SUBROUTINE NAME: 


HYBRID 


PURPOSE : 

This subroutine calculates transient temperatures using an optimum mix 
between implicit and explicit methods of solution. The explicit stability criteria 
of each diffusion node, CSG, is calculated on each temperature iteration as the 
capacitance divided by the sum of the conductors. This criteria is then checked 
against the user supplied time increment, DTIMEL. The temperature of these nodes 
with CSG less than DTIMEL are calculated using the implicit method of solution. For 
these nodes with CSG greater than DTIMEL, the explicit method of solution is used. 

T he order of calculations is arranged such that energy is conserved in 
conductors between the implicit and explicit nodes. Calculations are made on the explici 
nodes first. Next, the temperatures of implicit and arithmetic nodes are calculated 
using the latest explicit temperatures in the calculations. 

The implicit calculations are made using the methods described in Ref. 1. 
Using this method, temperatures of each node are made using the latest calculated 
adjacent temperatures. "Passes" are made repeatedly through the temperature calculation; 
until all temperature changes (between passes) have satisfied the user input tolerances 
DRLXCA and ARLXCA which must he supplied by the user. When the tolerance is satisfied 
for a node, the calculations of its temperature are temporarily suspended in the pass 
loop until all node tolerances are met. The calculation on all nodes are then resumed 
and the procedure is repeated until all node temperatures meet the tolerances on two 
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successive passes. The calculations may be over-relaxed or damped using user constants 
DAMPD and DAMPA. The default values for these variables are 1.0 for each. The maximum 
number of passes allowed through the temperature calculation loop is supplied by the 
user constant NL00P. Typical values for this variable are 500 to 1000, 

The implicit calculations for diffusion nodes may be backward difference, 
mid difference, or anywhere between backward and mid-difference. The first argument 
of HYBRID, ALPHA, determines the point in the iteration for evaluating the heat flux. 

If ALPHA = 1.0 or 0.0 (with a default value of 1.0) backward difference results. If 
0, < oi <_ 0.5, ALPHA is set equal to 0.5 and mid-difference results. If ALPHA is between 
0.5 and 1.0, the heat rate is ALPHA times that at the end of the iteration plus 
(1-ALPHA) times that at the start of the iteration. A second argument, K0P, will give 
a checkout print if f 0. Be prepared for a considerable amount of output if K0P f 0. 

The problem output is supplied at OUTPUT interval where OUTPUT is supplied 
as a user constant. The user may also supply a maximum allowable temperature change 
for the diffusion and ari thematic nodes by supplying values for DTMPCA and ATMPCA. 

If t u e changes are exceeded, the problem will be terminated. Default values for 
these are 1 . x 10®. 

RESTRICTIONS : 

The LPCS option is required and control constants TIMEND, OUTPUT, DTIMEL, 
NL00P, DRLXCA, and ARLXCA must be specified. Other control constants used or activated 
are: TIREN, TIME0, TIMEM, CSGMIN, DTIMEU, DTMPCA, DTMPCC, ATMPCA, ATMPCC, DAMPD, DAMPA, 
DRLXCC, ARLXCC, L00PCT, BACKUP, OPEITR, LINECT, PAGECT. 

CALLING SEQUENCE : HYBRID(ALPHA, KOP) 

DYNAMIC STORAGE REQUIREMENTS : 

This routine utilizes two dynamic storage core locations for each 
temperature node for non-flow problems or three dynamic storage locations for 
each temperature node for fluid flow problems. 
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SUBROUTINE NAME: 


INVRS 


See description for usage of SINVRS. 


SUBROUTINE NAME: 


QCOMB or ACOMB 


PURPOSE : 

QCOMB and ACOMB sum the interpolated value of the dependent variables 
of two arrays, A1 and A2, after multiplying A1 by ct-j and A2 by a 2 to form a third 
array, A3. For QCOMB, A3 contains all the independent variable values of both A1 

and A2 except where these values are equal. For ACOMB, the combined array will 
contain the independent variables of the A1 array only. 

RES T RICTIONS : 

Adequate space must be set aside in A3 but the space isn't required 
to be the exact amount needed by A3. 

CALLING SEQUENCE : QC0MB(A3,ctl ,A1 ,c*2,A2) or AC0M8(A3,al ,A1 ,a2,A2) 

where A3 is a doublet array with dependent variable values given by 

A3(i) = al*Al (i ) + a2*A2(1) 

al and aZ are constants to be multiplied times values of 
A1 and A2 at each point of A3 


Al and A2 are doublet arrays 



SUBROUTINE NAME: 


RADIR 


i 


PURPOSE : 

RADIR calculates the script-F values for infrared radiation heat transfer 
within an enclosure and uses these values to obtain the heat transfer 
during the problem. Several temperature nodes may be combined on a 
single surface for radiation heat transfer purposes. Also, the user may 
analyze problems with specular, diffuse or combinations of specular and 
diffuse radiation. See Section 2. 1.8.1 for definitions and detailed 
description of methods. 

RADIR calculates the script-F values on the initial call. This is per- 
formed by the procedure outlined in Section 2. 1.8.1. These values re- 
place the EFT values in the SC array for future use. The heat flux 
values are then calculated on all iterations by: 

(1) Calculating the temperature of each surface 

(2) Calculating the absorbed heat for each node 


The value given by equation 38 is added to the conductor sum for each 
node so that the proper convergence time increment may be obtained. As 
many enclosures as desired may be analyzed by each enclosure but each en- 
closure requires a different call to RADIR. RADIR must be called in VARIABLES 1. 

RESTRICTIONS : 

Must be called from VARIABLES 1 Surface nodes must be boundary nodes 
CALLING SEQUENCE : 

RADIR (A ( IC ) , SIGMA, TZERO) 

Where A is of the following format: 

A{IC),SN,SE,SR,SC,NA,SP,END 

SN,SE,SR,SC,NA, and SP are actual array numbers input using the *A procedure 
and are of the following formats 

SN(IC),n,SNl ,SA1 ,NN1 ,SN2,SA2,NN2 SNn ,SAn ,NNn,END 

SE(IC),SE1 ,SE2 SEn,END 

SR(IC) ,SR1 ,SR2 SRn,END 

SC(lC) ,SNF1 ,SNT1 ,EFT1 ,SNF2,SNT2,EFT2, — SNFm,SNTm,EFTm,END 
NA(IC) ,NN0(1 ,1 ) ,AN(1 ,1 ) ,NN0(1 ,2) ,AN(1 ,2) — NN0(1 ,NN1 ) ,AN(1 ,NN1 } 

NN0(2,1 ) ,AN(2,1 ) ,NN0(2,2) ,AN(2,2)~NN0(2 ,NN2) ,AN(2,NN2) 

r til t i 

i ill I l 

t ill l t 

NN0(n,l ) ,AN(n,l ) ,NN0(n,2) ,AN(n,2)— NN0(n,NNn) ,AN(n ,NNn) ,END 
SP(IC}, SPACE ,NSPACE, END 


The following definitions 
A 

SN 

SE 

SR 

SC 

NA 

SP 

n 

SN1 ,SN2, , . .SNn 

SA1 ,SA2, SAn 

NN1,NN2,....NNn 

SE1 ,SE2 SEn 

SRI ,SR2, , . . .SRn 
SNF1 a SNT1 a EFT! 

NN0{X,Y) 

AN(X,Y) 

NSPACE 

n 

SIGMA 

TZERO 


apply in the above calling sequence: 

Array i den itif cation for the array which identi- 
fies the other arrays containing the data 

Array number for the array containing surface 
numbers and areas 

Array number for the array containing the surface emis- 
sivities (may not be used in more than one call to RADIR) 

Array number for the array containing the sur- 
face reflectivities 

Array number for the array containing the sur- 
face connections data 

Array number for the array containing the tempera- 
ture node numbers and areas 

Array number for the array containing the space 
which is used for obtaining script FA values and 
for subsequent temperature calculations 

The number of surfaces 

Node number for surfaces - must be boundary nodes 

Total area for each surface 

Number of temperature nodes on each surface 

Emissivity values for each surface 

Diffuse reflectivity values for each surface 

Connections data: Surface number from, surface number 
to, E value from SNF1 to SNT1 , etc.(SNFl f SNT1 ) 

Temperature node numbers on surfaces; Node 

number Y on surface X 

Area of node Y on surface 1 

Number of spaces needed to «• tore script-FA 

values.- NSPACE must be an integer values 

of n *n(n+l)/2 

The number of surfaces 

Stef an- Boltzmann constant 

Temperature of absolute zero in problem units 
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SUBROUTINE NAME: 


RADSOL 


PURPOSE : 

RADSOL calculates a pseudo script-F for radiation from an external source 
entering an enclosure and uses these values to calculate the net heat 
transfer to each node due to the entering sou^e. A number of temperature 
nodes may be combined on a single surface for radiation purposes. Also, 
problems with specular, diffuse, or combinations of specular and diffuse 
radiation may be analyzed. Section 2. 1.8. 2 should be consulted for defini- 
tions and descriptions of methods. 

RADSOL calculates the pseudo script-F values on the initial call, as described 
in Section 2. 1.8. 2. The values are stored in the EFT values of the SC array 
supplied by the user. The heat flux values are then calculated on each 
iteration. 

The user may analyze as many enclosures as desired by supplying a call statement 
for each enclosure. Also, a user may analyze several wave length bands by 
supplying a call to RADSOL for each wave length band. 

RESTRICTIONS : 

Must be called From VARIABLES I; Surface nodes must be boundray nodes 
CALLING SEQUENCE : 

RADSOL (A( IC) ) 

Where the A array is of the following format: 

A( IC) ,SN,SE,SR,HT,SC,NA,SP,END 

SN,SE,SR,HT,SC,NA, and SP are actual array numbers input using the *A pro- 
cedure and are of the following formats: 

SN( IC) , n,SNl ,SA1 ,NN1 ,SN2,SA2,NN2 , -Snn ,SAn ,NNn »END 

SE(IC) ,SE1 ,SE2 , SEn,END 

SR(IC) ,SR1 > SR2,-- l SRn,END 

HT(IC) ,SHT1 ,SHT2 SHTn ,END 

SC(IC) ,SNF1 ,SNT1 ,EFT1 ,SNF2»SNT2,EFT2, — SNFm,SNTm,EFTm,END 
NA(IC) ,NN0(1 ,1 ) ,AN( 1 ,1 ) ,NN0(1 ,2),AN{1,2) — NN0(1 ,NN1 ) ,AN(1 ,NN1 ) , 

NN0(2,1 ) , AN (2,1 ) ,NN0(2 ,2) ,AN(2 ,2)— NN0(2,NN2) ,AN(2,NN2) , 

I I I I t t 

til* I I 

1111 I 1 

NNG{n ,1 ) ,AN(n»l ) ,NN0(n,2) ,AN(n,2) — NN0(n,NNn) ,AN(n,NNn) ,END 
SP(IC), SPACE, NSPACE, END 
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The following definitions apply in the above calling sequence 


A 

SN 

SE 

SR 

HT 

SC 

NA 

SP 

SN 1 ,SN2, , . .SNn 

SA1 ,SA2, . , .SAn 
NNl,NN2,...NNn 
SEl,SE2,...SEn 
SRI ,SR2 , . . . SRn 
SHT1 s SHT2, . . .SHTn 

SNF1 ,SNT1 ,EFT1 

NN0(X,Y) 

AN(X,Y) 

NSPACE 


n 


Array identification for the array which 
identifies the other arrays containing the 
data 

Array number for the array containing surface 
numbers and areas 

Array number for the array containing the surface emis- 
sivities (may not be used in more than one call to RADSOL) 

Array number for the array containing the 
surface reflectivities 

Array number for the array containing the 
incident heat curves or constant heat flux 
values 

Array number for the array containing the 
surface connections data 

Array number for the array containing the 
temperature node numbers and areas 

Array number for the array containing the 
space which is used for obtaining script values 
and for subsequent temperature calculations 

Node number for surfaces ; must be boundary 
nodes 


Total area for each surface 

Number of temperature nodes on each surface 

Emissivity values for each surface 

Diffuse reflectivity values for each surface 

Incident heat flow on surfaces; may identify 
curves containing incident values vs time 

Connections data: Surface number from surface number to, 
E value from SNF1 to SNT1 , etc. (must include SNF1=SNT1) 

Temperature node numbers on surfaces: Node 

number Y on surface X 

Area of node Y on surface X 

Number of spaces needed to store script- FA 

values - NSPACE must be an integer values of 

n(n+l)/2 

Number of surfaces 
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SUBROUTINE NAME : REVPOL 

PURPOSE : 

This subroutine performs single variable linear interpolation 
on a doublet array of X»Y pairs in the same manner as D1DEG1 except in re- 
verse order. The array is interpolated in reverse order to obtain the value 
of independent variable, X, which corresponds to the input dependent variable, 
Y. 


RESTRICTIONS : 

All values must be floating point numbers. 


CALLING SEQUENCE : 

REVPOL (Y,A(IC),X) 

where Y - input value of dependent variable 

A - Doublet array of X,Y pairs 

X - output value of independent variable 


SUBROUTINE NAME: 


SINVRS or INVRS 


J 


PURPOSE : 

These subroutines perform matrix inversion for symmetric, positive 
definite matrices using the efficient Square-Root or Symmetric Cholesky method. 

This method requires approximately half the computer time to obtain an inverse 
using the Gauss Elimination and Gauss-Jordan methods. Also, a significant increase 
in the accuracy has been observed. The symmetric matrix may be stored in half the 
square matrix space if desired. The Inverse is returned in its original space 
of the A-matrix. 

For SINVRS, the A matrix may be either a full square matrix or the 
upper triangular half of a square matrix. A check on the integer count relative 
to the matrix size is used to determine whether 1/2 matrix or full matrix is stored. 
INVRS assumes only the upper traingle of the symmetric matrix is stored. The (1,1) 
element is stored in the third data value of A for SINVRS and in the first data value 
of INVRS. The first data value of A contains the matrix size for SINVRS. 
RESTRICTIONS : 

The half syircnetric matrix must be stored as shown below for INVRS and 
for the half symnetric matrix option of SINVRS. Subroutine INVRS contains no error 
checks and should be used with extreme caution. 

CALLING SEQUENCE : SINVRS(A( IC) ,D) or INVRS(A(DV),N,D) 

where A is the matrix to be inverted and also, the inverse upon return 

D is the determinant of the original matrix to be inverted 

N is the matrix size 

The formats for A are as follows for SINVRS: 

(A) Full symetric matrix 

IC, N, BLANK, A(1 ,1) ,A{1 ,2) , - - - A(1,N) 

A{ 2,1 ) ,A(2,2) , A(2,N) 

t I I 

I I I 

I 1 t 

A(N,1),A(N,2) A(N,N) 
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SUBROUTINE NAME : TIMCHK 

PURPOSE : 

Subroutine TIMCHK compares the elapsed computer time against the 
requested computer time, RTIME, and terminates the run if RTIME is exceeded 
by the elapsed time. If the second argument, KODE, is non-zero an output 
of computer time used will be printed out on each call to TIMCHK. Thus, a 
call to TIMCHK in VARIABLES 2 should normally be with K0DE=0. If the output 
of computer time used is desired, TIMCHK should be called from OUTPUT with 
KODE f 0, The most desirable procedure is to supply two calls to TIMCHK : 
(!) a call in VARIABLES 2 with KODE = 0 and (2) a call in OUTPUT with 
KODE f 0. 

RESTRICTIONS : 

KODE should zero when called from VARTABLES 1 or 2. 

CALLING SEQUENCE : 

TIMCHK (RTIME, KODE) 

where RTIME = maximum computer time requested 

KODE 51 print code: = 0, computer time used is not printed out 

f 0, computer time used is printed out on 
each call to TIMCHK 
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SUBROUTINE NAME: WPRINT 

PURPOSE : 

Subroutine WPRINT will write all the values of the flowrates, 
pressure drops, pressures and valve positions. All values are printed out versus 
the actual numbers for which they occur. 

Should be called from OUTCAL 
WPRINT CK1,K2,K3,K4) 

K1 = 0, no flowrates will be printed 
= 1, flowrates will be printed 

K2 = 0, no pressure drops will be printed 
1, pressure drops will be printed 

K3 = 0, no pressures will be printed 
1, pressures will be printed 

K4 = 0, no valve positions will be printed 
1, valve positions will be printed 

DYNAMIC STORAGE REQUIREMENTS : 

Dynamic storage required by WPRINT is NW + NPR + NV where 
NW is the number of tubes, NPR is the number of pressure nodes, and NV is the 
number of valves. 


RESTRICTIONS : 

CALLING SEQUENCE : 
where: 


I 
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6.0 SAMPLE PROBLEM 

A sample problem was prepared for the SINFLO routine to demonstrate 
the input and output for a typical thermal/flow analysis problem. A schematic of 
the problem is shown in Figure 9 . The problem consists of 8 two dimensional 
radiator panels, each modeled by two flow paths (one for the main panel of 11 
tubes and one for the prime bypass tube). Contained in the system are a pump, 
a bypass valve (valve No. 1) and a stagnation valve between the two flow paths. 

The heat load to the radiator system comes through a counter flow heat exchanger which 
has a controlled inlet temperature of 40°F. The fluid is Freon 21 in the radiator 
system and water on the cooled side of the heat exchanger. The nodal subdivision 
for the fluid system is shown in Figure 9 . The structural nodal subdivision is 
shown in Figure 10. 

The sample problem was analyzed using the SNFRWD solution routine. 

The input for the problem is listed in Table 6 and the printed output is listed 
in Table 7. A few selected items were plotted using the plot package described 
in Appendix C. The plots of these items are presented in Figures 11 thru 17. The 
same sample problem was analyzed using the other temperature solution methods: 

CNFWBK, CNFWRD , CNFAST, HYBRID, CINDSS, SNFRWD FWDBCK, SNFRDL, CINDSL, and STDSTL. 
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3984,8 

Op 

16“ 

3909,8 

OP 

17“ 

42.884 

OP 

18. 

5 1 , 51 5 

Dp 

21“ 

337.79 

OP 

22* 

337. 79 

DP 

23- 

337.51 

DP 

24“ 

126.99 

Dp 

27* 

51.516 

OP 

28. 

**2,Ba3 

Op 

31- 

3980.6 

OP 

32* 

3980,6 

OP 

33. 

42,087 

DP 

~36- 

1 1288, 

OP 

’ 37- 

210.07 

OP 

38- 

737.75 

P 

1“ 

11708. 

P 

2* 

1 |499. 

p 

3- 

11289, 

P 

4* 

3)62.8 

P 

7* 

28)9,6 

p 

B. 

20U4, 0 

P 

11“ 

10839. 

P 

12“ 

6849,4 

p 

13- 

*797,9 

P 

1 4“ 

2775.6 

P 

17. 

2437,8 

p 

18. 

23)0,8 

P ' 

21“ 

6703.5 

P 

22* 

2722,9 

p 

23. 

2 1 0*07 

VP 

2. 

,33304 

VP 

3« 

,99999 





COMPUTER Time * .609 MINUTES 

The PROBLEM IDENTIFIED A 5 0»243“ has bEEr STORED AT ThIS POINT 
end of data 

• .Divide Check has occurred. • 


SINqa 


UN| V AC" 1 lOR FOFTHAN-v VERSION 


PAcE 


3 


13K’ - ’P-Pr 


T 

216* 

12.943 

T 

217“ 

7.9432 

T 

218- 

12*153 

T 

222* 

12.483 

T 

223* 

7.4592 

T 

224* 

11.903 

T 

226" 

-20.046 

T 

229“ 

-74,456 

T 

2 30* 

• 1 5 .856 

T 

234“ 

-20,120 

T 

235“ 

-24,463 

T 

236* 

-1&.971 

T 

2*)U“ 

-20,103 

T 

241* 

-74.447 

T 

242“ 

-15.952 

T 

246* 

-20,047 

1 

247“ 

-74^4*1 

T 

240“ 

-lS .893 

T 

252* 

64.B47 

T 

253“ 

64.577 

T 

254“ 

*7 .500 

7 

258* 

64.857 

T 

259“ 

64.538 

T 

2*0“ 

67,504 

T 

264- 

64,857 

T 

2*5“ 

64.638 

T 

266“ 

67,504 

T 

270" 

64.847 

T 

27 1“ 

64.577 

T 

272- 

67,501 

T 

276- 

63,579 

T 

277- 

63,2*4 

T 

278» 

6*. |97 

T 

282 - 

63,589 

T 

283“ 

*3.275 

T 

2®4- 

*6.200 

T 

288. 

*3.589 

T 

289» 

63^274 

T 

290* 

6*. |99 

T 

294- 

*3,579 

I 

295“ 

63.2*3 

T 

296“ 

6 * , 1 9* 

T 

30U* 

*8,827 

T 

3l)l* 

12.031 

T 

302“ 

12.032 

T 

306- 

— 1 5 , 9q 1 

I 

307“ 

*8 *825 

T 

308“ 

60,826 

T 

3U“ 

*7,503 

f 

3 1 3“ 

*6 1 1 98 

T 

31 4“ 

*6. |9B 

T 

400- 

-459,69 







A 

4. 

416.38 

a 

5“ 

209,35 




W 

9- 

416.26 

« 

lU“ 

83?, 64 




a 

14" 

416,15 

A 

15“ 

417.51 




w 

19" 

832.64 

A 

20“ 

416,42 




Vi 

24- 

206,93 

fl 

25- 

4)6,73 




W 

29- 

417.52 

A 

3 0 * 

416.15 




A 

34- 

16*7.3 

A 

35- 

25P0 .0 




* 

39. 

1*67.3 







Dp 

4* 

12,724 

Op 

5“ 

341,76 




DP 

?■ 

12.922 

DP 

lD“ 

I&U27 




OP 

14" 

3984.8 

OP 

15“ 

39)14.8 




DP 

19* 

15,527 

OP 

2d“ 

I2>31 




OP 

24* 

337.5) 

DP 

25- 

13.210 




OP 

29. 

3900,6 

OP 

3o“ 

39*0,6 




OP 

34 . 

412.15 

DP 

35“ 

21 no . 7 




DP 

39 . 

0154.2 







P 

4 . 

3175.5 

P 

5“ 

2832.5 




P 

9 . 

10877. 

P 

lu» 

6892.3 




P 

14" 

2708.5 

P 

15“ 

24c | .0 




P 

19“ 

6746.4 

P 

2fl" 

77*5.8 




P 

24“ 

,00000 
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SNFRBj^ PLOTS 


PLOT P S 0 0 R * I) 

title - sinflo sample probleh / siif1i«d 


• OOO HRS. TO 3*000 H«S» • PITH 3<00tl HRS'. pER g«IO 



the history Tape label is " 

SAHPl E PROBLEM / snfrao 

ThE JTEh courts are - 9AA 39&P 24pR jyP OBB UCC 

03T 39ffl 0F T otT 233ST 


ft 



1 

ITEM 

-9«ST 

*T 

320 

2 

ITEM 

I0*ST 

AT 

32S 

3 

iteh 

IHST 

AT 

33b 

<1 

ITEfl 

U55T 

AT 

'337 

S 

item 

-l I7ST 

AT 

339 

A 

item 

200ST 

AT 

342 

7 

item 

I96ST 

AT 

340 

a 

item 

199ST 

AT 

341 

9 

item 

-IFR 

AT 

67 

to 

item 

2FH 

AT 

60 

it 

item 

JbFfi 

AT 

102 

12 

item 

-2FR 

AT 

6U 

13 

item 

3FR 

AT 

69 

H 

item 

1 1 FR 

AT 

77 

IS 

item 

-|PR 

AT 

41 

lb 

item 

2PR 

AT 

42 

1 7 

item 

3PR 

AT 

43 

IS 

iteh 

24PR 

AT 

64 

1* 

ITEH 

*<ipn 

AT 

44 

20 

item 

9Pfl 

AT 

49 

21 

item 

IflPR 

AT 

5a 

22 

item 

2JPR 

AT 

63 

23 

item 

«2 VP 

AT 

65 

2T 

item 

3^P 

AT 

66 



I 



PLOT RUN PRINTED OUTPUT 


SNFBIt; PLOTS 


7 &.SS7 

.035 

23ni.732 

76.297 

.025 

2303.390 

74.953 

.025 

2305.295 

73.529 
.023 
23Q7. 140 

72. 189 
.□25 
23 n 6. 412 

70.979 

.025 

2309.744 

*9,909 
• 025 
23)0.700 

60.942 
,025 
23 i 2 • 33 7 

60,231 

,025 

23 t 3,335 

*7,722 

26.639 

2293.25? 

67.263 

60 .336 
2249.330 

47.101 

45,665 
2162.94 i 

*7,150 
I JO. 360 
2151.206 

47.117 

155.25* 


positioning *hd reading the history tape 


65,463 

70.00D 

2499,9*0 

2499,7)3 

209,330 

1 .000 

60,030 

70.000 

2499,9*7 

2499,710 

209,390 

1.000 

54,635 

70.000 

2499,969 

24’9,7l9 

209,556 

1 .000 

50.123 

70.000 

2499,960 

2499, 7 1 T 

209,727 

1 .uoO 

45,923 

70.000 

2499,970 

2499.719 

209.644 

1.000 

42,214 

70.D00 

2499,973 

2499,720 

209,972 

1 *000 

30,972 

*9,999 

2499,973 

2499. 7| 9 

2 1 U .057 

1 .000 

36, | 10 

69,999 

2499,970 

2490.731 

210.201 

,900 

33.590 

69,539 

2499 ,969 

2373.3*9 

210.296 

.949 

31,330 

*7,340 

2473,356 

22*5.130 

2I2.Q35 

,9|3 

29,295 

*9,209 

2439.650 

2136.312 

212.50* 

.072 

27.407 

it .093 

2404.329 

1994,45* 

211.057 

.028 

25,607 

60,449 

2371,620 

1065.030 

212.440 

.763 

23. R** 

59,002 

2344,740 

1730.019 


LOflnEo 

.01030 

66.1*7 

* 6 * 5 1 5 

. 2»0 

11196.5*6 

1 ,OuO 

LOAOtO 

.02000 

*0.963 

61.522 

. 2*0 

U 179.351 

1 .OUQ 

LOADtO 

,03000 

53,74) 

56. 360 

4 2*0 

1 1 144.036 

1 .OuO 

LOADtO 

.04000 

51 .QuS 

51.553 

.247 

11131,114 

1 , Qua 

LOADtO 

.05000 

46.715 

47,211 

• 250 

11 1)9.096 

1 « 0 uO 

LOAOtO 

. 0*000 

42.917 

43.356 

,250 

11119,116 

1 , 0 uO 

LOAOtO 

.07000 

39.593 

39 1 90 | 

,240 

I 1 100.500 

1 , OuO 

LOAOtO 

.00000 

36,6/1 

37.022 

5 l , 2 J9 

I 1099,77; 

1 .OLD 

LOADtO 

,09ooo 

34,721 

34.972 

126,605 

107*9,003 

1 , 0 uO 

LOAOtO 

•lOOOQ 

33,4*2 

33.621 

2l4. 2 U 1 

10293.401 

,969 

LOAOtO 

. 1 I 00 Q 

32,SaS 

33.054 

313,4*7 

9646.96* 

.976 

LOADtO 

. 1 2000 

32.047 

32.924 

413.05) 

052U.132 

.9*2 

LOAOtO 

• 1 3000 

31.7*1 

33.134 

515. * u 7 

0000.927 

.949 

LOAOtO 

•1400U 

31,742 

33.599 

614,792 

7327.946 


HRS,' LOOKING for 

7 7 , 23q 60.000 

10970.103 lo7*3.l99 

hRS. LOOKING FOR 

76^2*3 0O.DOO 

10»72,363 | q 742 • 6fl3 

HHS. LOOKING for 

^4 • 000 80.000 

10935,638 1Q726.4SI 

HRS, _ LOOKING FOR 
73 jRS? 0O.UDO 
10921.897 Iq712,*73 

HRS. LOOKING FOR 

72^11. Bo.uoo 

10909,845 ti.70ll.59S 

HRS, LOOKING FOR 

70^914 60.000 

10909,0*2 lo*90. 007 

HRS, LOOKING FOR 

69,046 00.000 

IO09|J. 47q 10601,174 

CHS. LOOKING FOR 

46.91! 00.000 

10713,397 10350.370 

Hfl5, looking FOR 

tfi;i9 z ao.ooo 

1030*. 660 9076,909 

HRS, LOOKING FOR 

67 j 49<t 00,000 

9700,077 9220.994 

HRS, LOOKING FOR 

67.259 00,000 

9097,71* 0450.561 

MRS, LOOKING FOR 

67 j09 1 0O.UDO 

7061.721 7590.742 

HNS. LOOKING FOR 

4711&J 00,000 

7496,114 *929,097 

HNS, LOOKING FOR 

67,115 00.000 

6021,377 4275.360 


3 .OUuOO 

HR 5 i 



7 1 .uao 

2*99,94* 

2499,9*6 


.OUO 

9924, 6*4 

2302.739 


3. UOuOO 

HHS, 



*7, <364 

2499,99* 

7499,9*7 


*000 

990*. 094 

7303,400 


3.00000 

HRS. 



63.370 

2499,997 

7499,9*9 


,000 

9669,041 

7305,313 


3.00uU0 

HHS. 



SB. 731 

2*99,997 

7499,968 


.000 

9675,9*2 

7307.159 


3.00U00 

HRS « 



54,397 

2500, 0(iU 

7449.970 

_3 

• OUO 

9663.773 

73Q0.632 

> 




a 

3.00000 

hR5 ^ 


r 

50.440 

2*99,995 

7499,973 

pi 

• UuO 

9 B 53 ,0V U 

7309,706 

CD 

3 • OOuUO 

hRS . 



47,0 22 

2499,996 

7499 . 973 

O 

.OUQ 

9B44.191 

7310.001 

u 

Z 




*-j 

3. OOUUO 

HRS. 


5 

43,403 

2499.997 

7499,970 

a 

• COG 

9442,904 

7322.517 

n 

3 .OOuUO 

hRS. 



41,319 

2499,99* 

7499,969 


,uu(l 

9 112.537 

7372.023 


3, OOUUO 
39,475 

HHS. 

2499,99* 

7473.356 

0 0 

»uUO 

6422.790 

7457.562 

3.00000 

HRS. 


30. <59 

2499,997 

7439,656 


.000 

7612,693 

7607, O 5O 


3, OOiJLrO 

hNS. 



37.749 

2499,997 

7404,329 

.000 

7n 34.300 

7757,170 

3.0UUU0 

hHS. 


S 

37. *3 t 

2*99,990 

7171 ,628 

3$ 

< UUO 

* *i 1 * . b 7 j 

TO jS ,394 

3, OOUUO 

HflS, 



37, *20 

24*9,999 

7344,740 


.OUO 

563 | ,U4 „ 

1350,738 




SHFR« d PLOTS 

43*825 -l&. 901 44.1*0 J9.4l7 3**419 *025 80,000 *0.0*2 2«99.9*/ ?h95,9a* 

♦ 026 2H**,*6* 832*493 1447.323 117q7,7*o HHVS.SOn 11299.206 *uOO 3)76. 5j*l |f>8T7.U!> 

23 | 0*798 210*076 .3)3 l»Oufl 

25 OAtA VALUES HAVE 0 EEn STORED FOR EACH OF |ol T J HE POINTS 

2500 data values have seen stoned 



TABLE a (CONTINUED) 


SNFRWp PLOTS 


...... 



ITEM 

TYPE 

a vg plotting symbol and description 



T-HIN 

Y-HAK 

STATUS 

1 

-9a 

ST 

l 

RADIATOR INLET TEhPCNATORE 


W ® 

OEGF 

-1 *590*01 

7 <686.01 

320 

2 

106 

ST 

2 

MAIN RADIATOR OUTlET TEMPERATURE 

- - 

oegf 

u.Uuu 

0.000 

328 

3 

119 

' ST 

3 

phihe tube Outlet temperature 

m m 

OEGF 

0 tQOO 

0.000 

336 

9 

IIS 

ST 

8 

HIKED RAOTATQH OUTLET TEMPERATURE 


OEGF 

U.Ouu 

0.000 

337 

5 

-117 

ST 

1 

RADIATOR CONTROLLED OUTLET • 

wx inlet 


OEGF 

3. 292*m 

8.00Q*01 

339 

6 

200 

ST 

2 

HA OUTLET ON RAOIaTOh SIDE 


»• 

OEGF 

O.OoD 

U .000 

392 

7 

198 

ST 

3 

11 K INleT on RATER 5 1 qE 



OEGF 

o.Ooo 

0,000 

3M0 

9 

1 99 

ST 

9 

ha outlet on rater Side 


• a 

DEGF 

O.OGO 

0.000 

391 

9 

-1 

F R 

1 

total pump flow rate 

lb/hr 



2 * 3u?-r>2 

2.900*03 

67 

10 

2 

Fft 

2 

Total radiator flc« rate — 

LD/Hfi 



D * Qud 

U.000 

68 

. 1 

36 

F ft 

3 

BYPASS FLOW RATE 

lb/hr 



d « OuO 

0,000 

102 

1 ? 

-2 

FB 

i 

total radiator flow RATE — 

lb/hr 



2.971-01 

2,500*03 

60 

13 

3 

FB 

2 

MAIN RADIATOR FLOft RaTE — 

lb/hr 



OoOQU 

0.000 

69 

19 

! 1 

FR 

3 

prime tube Flow Rate 

LB/HR 



O.OUO 

0.000 

77 

IS 

-1 

PR 

1 

PUMP OUTLET pressure 

-- PSF 



u.ono 

1.171*09 

•M 

ia 

2 

PS- 

2 

valve 1 INLET PRESSURE 

— PSF 



0 i OuD 

O.OOO 

92 

17 

3 

PS 

3 

VALVE 2 INLET PRE&SUrL 

— PSF 



o.oou 

0 ,000 

93 

18 

29 

PR 

9 

PUMP INLET pressure 

— PSF 



o .Ona 

0.000 

69 

19 

-9 

Pft 

1 

main maDIATOR INleT PRESSURE— psf 



2.003*112 

l .060*09 

99 

20 

9 

PR 

2 

prime tube inlet pressure 

— PSF 



u .odd 

O.Qbd 

99 

21 

18 

PR 

3 

PRESSURE AT RADIATOR OUTLET 

— PSF 



d.Ood 

0.000 

SB 

22 

23 

PR 

9 

PRESSURE AT ha INlET 

— PSF 



d * Qou 

0.000 

63 

23 

-2 

VP 

1 

VALVE 2 POSITION 




3. 33 i -n l 

1 e DQO + OO 

65 

29 

3 

VP 

2 

VALVE 3 POSITION 




0 . 0 1. U 

0.000 

66 
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8 S 




0 ■ i 


O'.. 


0 I 

a J 


a 1 s 


SNFRlU PLOTS 

STARTING PLOTS 


PLOTTING 

1 

RADIATOR INLET TEHPERATuHE 


« m 

oer.F 

plotting 

2 

main radiator outlet temperature 

m m 

oegf 

plotting 

3 

PRIME TUBE OUTLET TEMPERATURE 

m • 

OEGF 

— plotting 

H 

MIXED RADIATOR OUTLET TtHPERATUflE 

•- 

OEGF 

PLOTTING 

1 

RADIATOR CONTROLLED outlet. 

HX inlet 

““ 

OEGF 

plotting 

2 

HX OUTLET on RADIATOR side 


m m 

OEGF 

plotting 

3 

MX INLET ON HATER SICE 


w * 

oegf 

PLOTT [NG 

N 

HX OUTLE 1 ON WATER SIDE 


— 

OEGF 

plotting 

I 

total pump flop rate 

lb/hr 



— PLOTTING 

2 

TOTAL RAutATOR FLO" R*T E — 

LB'rtR 



PLOTTING 

3 

BYPASS FcO* R*TE 

lB/mR 



PLOTTING 

1 

TOTAL RADIATOR FLOH rATe — 

LO/hR 



PLOTTING 

2 

MAIN RADIATOR FlON R*TE — 

LB/Hfi 



PLOTTING 

3 

PRIME TUoE FLOW RATE 

L0/HH 



PLOTTING 

| 

PUMP OUTLET pressure 

— PS," 



PLOTTING 

2 

valve l inlet Pressure 

-- PSF 



flatting 

3 

valve 2 INLET PRESSURE 

-- P5F 



plotting 

B 

PUMP INLtT pressure 

— PSF 



PLOTTING 

t 

MAIN RADIATOR INLET PRESSURE*- PSF 



PLOTTING 

2 

PRIME TUDE INLET PRESSURE 

• * PSF 



plotting 

3 

pressure AT RADIATOR outlet 

— PSF 



"" ■ ■ PLOTTING 

R 

pressure at HX INLET 

— PSF 



PLOTTING 

1 

valv: 2 position 




PLOTTING 

2 

VALV. 3 POSITION 





COMPUTER time e olqHJS minutes 


ePHOt; 


BM5G>H end OR213E 


BORkPt print* 
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XABI.E 8 (CONTINUED) 


TEMPERATURE 


FIGURE 11 

RADIATOR TEMPERATURE PLOTS 



SINFLO SAMPLE PROBLEM /• SNFRWD 
RADIATOR INLET TEMPERATURE 
MAIN RADIATOR CUTLET TEMPERATURE 
PRIME TUBE CUTLET TEMPERATURE 
MIXED RADIATOR OUTLET TEMPERATURE 


(HOURS) 





















TIME - (HOURS) 
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FIGURE 14 

RADIATOR FLOW RATE PLOTS 


Si»;"LO r .C~LEt! 
TOTAL RADIATOR FLOW RATE -■ 
MAIM RADIATOR FLOW RATE - 
PRIME TUBE FLGW RATE 


/ rtlFRWD 
LB /HR 
LB/HR 
L8/HR 










PRESSURE 




1 


!_ 




FIGURE 15 

SYSTEM PRESSURE PLOTS 


SIN c l.P S.V^I.F PRFR'.FM / SNFP'.'D 
PUMP OUTLET ' JRE -- PSr 

VALVE 1 INLET IRESSURE -- PSF 

VALVE 2 INLET PRESSURE — PSF 

PUMP INLET PRESSURE — PSF 

TTTT7T fl TTTlT!*Tl*T- 




12000 


10000 


•000 


MOO 


*000 


2000 


TIME - (HOURS) 







PRESSURE 


FIGURE 16 


RADIATOR PRESSURE PLOTS 


SINFLO SAMPLE PRTPLEM / SNFRS'D 
[1] MAIN RADI A TOR iNLc T PRESSURE-- PS? 

[21 PP! •• T'T INLET •- — FT 

i 3 1 PRESSURE AT RaOIAIuR OUTLET -- PSF 

[41 PRESSURE AT HX INLET — PSF 



TIME - [HOURS) 


Original pagb ig 
poor quaiot 
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RADIATION INTERCHANGE ANALYSIS 


Capabilities have been incorporated into subroutines for use with SINDA to 
facilitate the analysis of radiation heat transfer in an enclosure. The capabilities 
include the ability to: 

(1) Analyze diffuse and/or specular Infrared radiation in 
an enclosure 

(2) Analyze diffuse and/or specular radiation from an external 
source for as many wave bands as desired 

(3) Consolidate several temperature nodes into a single surface 
to improve computational efficiency 

A radiation surface is defined as a group of temperature nodes which may be 
assumed to have identical radiating properties, angle factors and interchange factors. 

The subroutines account for the net radiation heat transfer between a number 
of surfaces due to the emitted radiation from each surface, reflected radiation from 
each surface, and radiation from any number of incident sources. The reflection of the 
energy originally emitted by another surface or from an external source may be either 
diffuse, specular, or any combination of the two. 


Emitted Radiation In A Cavity 

The radiosity of a surface is defined as the flux of infrared radiation 
leaving that surface with a diffuse distribution (according to Lambert's Law). That 
energy leaving a surface which has been reflected in a specular manner does not contri- 
bute to the radiosity of that surface. The incident infrared radiosity is denoted by 
the symbol H. The reflectance (1 - e) of a surface is separated into two components, 
the diffuse reflectance (p), and the specular reflectance (p s ). Here e is the emit-' 
tance of the surface and Is equivalent to the absorptance for long wavelength radiation 
With the angle factors (Fij) defined in the normal way, there exist similar angle 
factors which relate the geometrical ability of surface i to radiate to surface j by 
means of a mirror-like reflection from specular surface k. Reference to Figure A-l In- 
dicates the method of imagery which will enable the calculation of these reflected anglf 
factors. Here the angle factor to surface d is identical with the angle factor to 
the image of surface j. Also the angle factor is limited by the ability of surface i 
to "see" through the "window" of surface k. With the specular surface angle factors 
so defined, an interchange factor Eij is defined similarly to reference 8 as follows: 

-becedins page blank not mam 

A-l 



FIGURE A-l ILLUSTRATION OF METHOD USED TO DETERMINE SPECULAR SURFACE REFLECTED VIEW FACTORS 
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Here is the angle factor from i to j as seen in the specular surface 

k, ) is the angle factor from i to j as seen in the double specular 

reflection from k and There are an infinite number of possible combinations 

of these multi-ref lections. It is evident that the interchange factors account for 
the specularly reflected radiant flux from the reflecting surface. This 
portion of total leaving flux is not a component of the radiosity of that 
surface. The radiosity may be written 
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Now the interchange factors obey the reciprocity relation 
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Substitution into the equation for B results in 
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This equation represents a set of linear, simultaneous, inhomogeneous 
algebraic equations for the unknowns (Bj). The ..symbol 6-fj is the Kronecker 
delta function which is 1 when i -’j and is 0 when i f j. 

Note that the coefficients of Bj in equation fA-3) do not form a 
symetri c coefficient matrix since the off diagonal terms contain - 
This equation can be made symetri c by multiplying each equation by A-j/p-j. 
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This gives 
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Written in matrix form this equation is 
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Where E is a symetric coefficient matrix. 
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The solution is 
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The net heat transfer rate asorhed by surface i is given by 

q i = V t [ H, - ^ T. 4 ] 

Where H-f is given from equation (A~2) as 

H, = -1- C B-j ~ ‘“-T, 4 ] 


Substituting in for H-j gives 
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Substituting in for B-j from aquation (A-5) into aquation (A~6) gives 
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Since, in steady state, Qi - o , and T^ 4 - T- 4 for all i and j we can 
conclude that ns J 
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Making the above substitution in equation (A'^7) -gives 
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This equation gives the heat flux between surfaces. However, each surface 
can contain several nodes. The heat absorbed by for each node is determined 
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Where n - the node number on surface i 


Prior to each iteration, the temperature of the surfaces are determined by 
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Where nn - the number of nodes on surface i 


Since the heat transfer rate given by equation (A-9T depends on the node 
temperature, stability considerations must be taken into account. This is 
handled by storing the following relation into the array containing the 
sum of the conductors used for time increment calculation 
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Subroutine RADI R makes the calculations necessary to obtain Qn 
given by equation (A-9) and CO Nn given by equation (A-1T). The following is 
a summary of the calculations: 

A. The following are performed the first time through RADIR: 

1. From the user input values of E-jj, A-j , and , the E 
matrix given by equation (A-4) is formed. Only half 
of the symetric matrix is stored to save space. 

2. The E matrix is inverted in its own space to get 
with elements e-j 

3. The^ A* , values are determined from equation (A-8) 
and stored in the surface connections data. 

B. The following calculations are performed on each temperature 

iterations: 

1 . The temperature of each surface is calculated by equation 
(A-10). 

2. The heat absorbed for each node is determined using 
equation (A-9) and is added to the Q array. 

■ A-6 . ' 


The routine utilizes data used for obtaining ^ A-j j in step A as working 
space for step B, thus, maximizing space utilization. 

Radiation From External Source 

As with the internally generated radiation, the solar (or any other 
external source radiation) interchange factor is defined by 
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Where p* s is the solar specular reflectance of surface K 

F..(K) is the angle factor from i to j as seen in the 
specular surface k 

F..(K,/) is the angle factor from i to j as seen in a double 
specular reflection from j to / to k back to i 

The interchange factors as defined above accounts for the specularly flux 
reflected from the surface. Thus, since the specular component of the flux 
is assumed to go directly from surface i to surface j by the interchange 
factor, Ei j , this portion of the total flux is not a component of the radiosity 
for the i nternimedi ate surfaces (k and \ above). The radiosity of surface 
i is given by 
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is the incident energy 
)* is the diffuse reflectance 
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The energy incident upon a surface is given by 
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Where Sj is the energy directly incident on surface i from 
an external source 


Substituting equation (A-T2) into CA-13l ? multiplying by and simplifying 
gives the following relation for the radiosity 
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This set of n equations can be written in matrix form as 


E* B* = S 


Note that the equations are written so that E* is a symetric matrix, which 
has the solution for’ B* • 
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Where is the ijth element of the inverse of the E* matrix 
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The heat flux absorbed by the i th surface is given by 
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Combining equations (A~16), (A-17), and (A-18) gives 
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Then the absorbed heat flux is given by 
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Equation (A-21) gives the heat absorbed by each surface. However, each surface 
may contain several temperature nodes. The absorbed heat for each node is 
given by: 

Q* | q* (A-22) 

Where A is the area of the node 
n 


Subroutine RADS0L was written to make necessary calculations to 
obtain Qn given by equation (A-22). The foil owing is a summary of the 
calculations: 


A. 


B. 


The following calculations are made the first time through 
RADSOL: 


1 . 


2 . 


3. 


From the user input values of E • p. , and Ai , the 
E* matrix given by equation (A-15) is formed. Only one 
half is stored since E* is symetric. 

The E* matrix is inverted in its own space to get E*' 
with elements, e.! 


-1 


The 0 ^*. Ai values are determined from equation (A-20) and 
stored in the surface connections data. 

The following calculations are performed on each temperature 
iteration: 

1. The heat flux absorbed by each node is calculated by 

% = -ib £ 
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2. The net heat absorbed by this wavelength radiation is 
calculated for each temperature node on each surface by 
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This quantity of absorbed heat is added to the Q array for 
node n. 


Note that the user may specify subroutine RADSOL for as many bands of 
radiation from an external source as desired. A single call is required 
for each band. 


APPENDIX B 
FLOW DATA STORAGE 


The flow data which is input in the FLOW DATA block described in 
Section 4.1 is stored by the preprocessor in Tabled common arrays. These 
arrays will be included in the main processor phase routine and the routines 
generated from the four operation blocks (EXECTN 9 VARBLT, VARBL2, and OUTCAL) 
The arrays will be dimensioned in the main processor routine. The following 
is a list of the arrays: 

1. Flow Data : /FL0DAT/FL0W(ND) , where ND is the amount of space 

required for the flow data array. This array includes 
tube connections data and tube data for all 
systems, specified pressure nodes, valve data, pump 
data and enthalpy curve 

2. System Data : /SYSDAT/SYSTEM(1 5, NS) , where NS is the number of 

systems. Systems data include property data, solution 
parameters and specified pressures. 

3. Fluid Lump Type Data : /TYPDAT/TYPE( 10, NTP) , where NTP is the 

number of fluid types. 

4. Flowrates : /WD0T/W(LT), where LT is the largest input tube 

number . 

5. Pressures : /PRESS/P{LP) , where LP is the largest input pressure 

node number. 

6. Flow Conductors : /FLOWG/GF(LT ) 

7. Valve Positions : /VALVP/VP{LV) , where LV is the largest input 

valve number. 

8. Imposed Flowrates : /WDOTI/WI (LP ) 

9. Added Flow Resistances : /FLOWR/AFR(LT) 

TO. Pressure Drops : /DELTAP/DP(LT) 

11. Dimensions : /FDIMNS/NTYPE, NSYS, NTB, NP, NV, NFD 

where NTYPE is given the value NTP (above) 

NSYS is given the value NS 

NTB is given the value LT 

NP is given the value LP 

NV is given the value LV 

NFD is given the value ND 


The items to be stored in the above are discussed in more detail 

below. 

1.0 FLOW Array - 

The FLOW array contains flow data that is not easily addressable by 
the user. I tons contained are (1) the network and subnetwork connections 
data, (2) the tube data which includes the fluid lump/tube lump pairs and the 
fluid lump type, (3) the specified pressure nodes, (4) the valve data, (5) 
the network valve data, (6) the pump data, and (7) valve data locations. 

The format for storing each item is discussed below* 


1.1 Network and Subnetwork Connections 

The network and subnetwork connections data is stored in the following 
order for each network or subnetwork: 

I Cl , 'NAME! ’ , LOCPRT , LOCVl , NTB11 , NFRM11 , NT011, L0CD11 

NTB12, NFRM12, NT012, L0CD12 

i -i i i 

i i i i 

■ i t i 

NTBln, NFRMln, NTOln, LOCDln 
IC2, 'NAMEn' , L0CPR2, L0CV2, NTB21 , NFRM21 , NT021, L0CD21 

i i t i 

« i i i 

NTB2n , NFRM2n, NT02n, L0CD2n 

■ t i i 

i i i i 

it * i 

ICn, 'NAMEn ' , LOCPRn, LOCVn, NTBnl , NFRMnl , NTOnl , LOCDnl 

i t i ' 

i i t • 

NTBnn, NFRMnn, NTOnn, LOCDnn 
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where ICf is the integer count of the number of spaces in the connections 
data for the ith network or subnetwork 
NAMEi is the 4 character name of the ith network or subnetwork input 
on the heading card 

LOCPRi is the location in the flow data array of the specified 
pressure nodes for the ith network or subnetwork 
LOCVi is the location of the ith network or subnetwork valve data 
(which is an array of locations of the actual valve data) 

NTBij is the tube number of the jth tube of the ith network of 
subnetwork 

NFRMij is the "from 11 pressure node for the jth tube of the ith net- 
work or subnetwork 

NTOij is the "to" pressure node for the jth tube of the ith network 
or subnetwork 

LOCDij is the location of tube data (or subnetwork connections) for the jth 
tube of the ith network or subnetwork 

If LOCD > 0 it is the location of the tube data (fluid/ tube 
lump pairs and type no's.) 

If LOCD < 0 it is the location of the subnetwork connections 
data for tube j 

If LOCD = 0, the user is supplying the flow resistance for 
tube j in the added flow resistance array, APR 

A sort is required on the connections data for each network or subnetwork. 

The connections must be arranged so that for each pressure node, all NTO 
references for that node must occur in the list prior to any NFRM references. 

The four data values (NTB, NFRM, NTO, LOCD) must remain intact as a group 
during the sort. Tubes whose "from" node, NFRM, is not referenced as a "to" 
node, NTO, should come first in the connections data. 



T.2 Tube Data 

The tube data portion of the FLOW array contains the fluid lumps, 
fluid lump types and tube lumps for each tube. This data is referenced by 
the LOCD values in the connections data for each tube described in Section 
4.1. The format for the tube data portion of the FLOW array is: 

ICl, NFLMP1 1 , NTYPE1 1 , NTBLMP11 , — , NFLMPln, NTYPEln, NTBLMP1 n . 

i • i i i • i 

i i i i i i i 

i i < i i < 1 

ICn, NFLMPnl , NTYPEnl, NTBLMPnl , — NFLMPnn, NTYPEnn, NTBLMPnn 

where ICi is the integer count for the tube data for tube i (must be 
a multiple of 3) 

NFLMPij is the relative fluid lump number of the jth fluid lump 
in tube i 

NTYPEij is the type number of the jth fluid lump in tube i 
NTBLMPij is the relative tube lump number for the jth fluid 
lump in tube i 

Notice that NFLMP and NTBLMP are relative lump numbers. Thus, during storage 
these numbers must be converted from actual numbers which are input to 
relative numbers. 

1.3 Specified Pressure Node Data 

The specified pressure node data is a list of the pressure nodes 
whose pressures are not calculated. One such list exists for each network in 
the problem and may also exist for any subnetwork if it contains any specified 
pressures. The format for each specified pressure node list is: 

IC, NSP-i , NSP 2 - - - NSP IC 

where IC is the integer count which is also the number of specified 
pressure nodes in the network or subnetwork 
NSP-j is the ith specified pressure node 

T.4 Valve Data 

The valve data described in Section 4.3 is stored in the FLOW array 
The format for this, is slightly different for the different types of valves. 
For the rate limited valve the format is: 






IC, NV, NTS! , NTS 2, MODE, XMIN, XMAX, E, TSEN1 , TSEN2, DB, RF, RL 
For the polynomial valve it is 

IC, NV, NTS1, NTS 2,. MODE, XMIN, XMAX, E, TSEN1 , TSEN2, CO, Cl, CH, C3, C4, C5, VTC 
For the switching valve it is 

IC, NV, NTS!, NTS 2, MODE, XMIN, XMAX, E, NSEN, T1 , T2 

where the symbols are described in Section 4.3 The integer count for each 
is the number of data values and is 12 for a rate limited valve, 16 for the 
polynomial and 10 for the switching valve. 

1.5 Network Valve Locations 

The network valve locations is a list of locations in the FLOW array 
for the valve data of the valves in a network or subnetwork. One such list 
is needed in the FLOW array for each network or subnetwork that contains valves. 

The location of the network valve data list is provided in the fourth location 
of the Network Connections Data. 

The format of the network valve locations is: 

IC, L0CV1, LOCVIC 

where IC is the number of valves in the network 

LOCVi is the location in the FLOW array of the valve data for 
the ith valve in the network 

1.6 Flow Source Data 

The flow source data lists are supplied in the FLOW array for each 
flow specification statement input in the BCD FLOW SOURCE data block described 
in Section 4..4.. The input statements from the FLOW SOURCE data are transfered 
directly to the FLOW array except an integer count is added to each list and 
array numbers are converted from actual to relative numbers. 

The formats for storage are as follows for the three types of flow 

sources: 

Flow As A Function of Time 
IC, NPI, AW 

Pressure Rise As a Tabulated Function of Flowrate 
IC, NPI, NPO, ADP 

Pressure Rise As A Polynomial Function of Flowrate 


IC, NPI, NPO, CO, Cl, C2, C3 , C4 


where IC is the integer count of the list (2, 3, and 7 respectively) 

All other variables are described in Section 4.4, The actual numbers of arrays 
referenced by AW and ADP must be converted to relative locations prior to 
storage in the FLOW array. Only one flow source data list per network is to 
be stored in the FLOW array and the location is referenced from the SYSTEM 
array (to be discussed later). 

An option on AW is that it may be input as an array or a real con- 
stant. If AW is supplied as a real constant, the flow source list is not 
stored in the FLOW array. Rather, the constant, AW, is stored in the imposed 
flowrate array, WI, (to be discussed later). 


1 , 

1.7 Valve Locations 

The valve locations list is a list of locations of the valve data 
(whose input is described in Section 4.2 and storage is described in the 
Appendix) for all the valves in the problem in order of valve number. There is 
only one valve location list in the FLOW array and the location of this list 
is given as the seventh item in the FDIMNS Tabled common block (described be- 
low). 

The format for storage of the valve locations is 
IC, L0CV1 , L0CV2, , LOCVIC 

where IC is the interger count and is the total number of valves in the 
probl em 

LOCVi is the location in the FLOW array for the valve data for valve 
number i 

2.0 SYSTEM Array 

The system array contains fluid property data (or locations of 
property data), the gravitational constant (gc), solution parameters, and 
the locations in the FLOW array for the flow source list, the network connections 
data and the enthalpy curve for each system. The SYSTEM array is a two 
dimensional array dimensioned to 15 by NSYS where NSYS is the number of systems. 
Thus 15 locations are allocated for each system (only 13 are currently used 
leaving 2 blank spaces per system). The system array is in the Tabled common 
block SYSDAT. 
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The format for storage of the SYSTEM array is 




ACP-j jARO-j ,AMU 1 iAKT-j ,8^ .MPASS-, ,T0L 1 .MXPASS-j ,FRDF-| ,K0P-j ,LOCP 1 .LOCNET, ,L0CH 1 ,0,0 

acp 2 ,aro 2 ,amu 2 ,akt 2 ,gc 2 ,mpass 2 ,tol 2 ,mxpass 2 ,frdf 2 ,kop 2 ,locp 2 ,locnet 2 ,loch 2 ,o,o 


I I I 


ACP n* ARO n sAHU n’ AKT n ,GC n ,HPASS n ,TOL n ,MXPASS n ,FRDF n ,KOP n ,LOCP n yLOCNET n ,LOCH n i0,D 


where ACP.. ,AR0.. ,AMIL , and AKT. are the relative array numbers for the 

arrays or the values of the constant values 
for the specific heat, density, viscosity and 
thermal conductivity for the ith system 
GC is the gravitational constant for the ith 
system 

MPASS. is the number of temperature iterations be- 
tween pressure solutions for system i 
MXPASS^ is the maximum number of passes in the 

balancing loop permitted to obtain a pressure/ 
flow solution for system i 
FRDF.J is the flowrate damping factor for system i 

T0L.J is the solution tolerance on the fraction 

of change of flowrates from one pass in the 
flow solution to the next for system i 
KOP. is the check-out-print code for system i 

LOCP. is the location of the flow source data in 

the FLOW array for the ith system 
LOCNET^ is the location of the network connections 
data in the FLOW array for the ith system 
L0C1-K is the location of the enthalpy curve in 
the FLOW array for the ith system 


The values for ACP il ARO.,AMU.,AKT i ,GC.,MPASS i# MXPASS i ,FRDF i JOL. and KOP. are 
taken from the systems input supplied in the BDC 3NETW0RK block except that 
array numbers are converted to relative array locations for ACP, ARO, AMU 
and AKT and default values are supplied for GC, MPASS, MXPASS, FRDF, TOL and 
KOP if no values are input {Default values are shown in Section 4 t ll. The 

values for LOCP, LOCNET and LOCH which are storage locations in the FLOW array 
are determined as the FLOW array is built during the preprocessor phase. 



3.0 


Fluid Type Array 

The fluid lump type data is stored in the TYPE array which is in 
the TYPDAT labled common block. This array contains the fluid lump type 
information which is input in the BCD 3FLUID LUMP DATA input block on the left 
of the equal sign for all type cards. The TYPE array is a two dimensional 
array, dimensioned to 10 by NTP, where NTP is the number of types. The format 
for the TYPE array is 

CSA-j ,WP 1 ,FLL 1 ,AHT ] ,NHL 1 ,MFF 1 ,FFC ] s Fi^ , FLL^WP., 

4.0*CSAf 

csa 2 ,wp 2 ,fll 2 ,aht 2 ,nhl 2 ,mff 2 ,ffc 2 ,fi 2 ,f2 2 , fll 2 *wp 2 

4.0*CSA ? 

i i i i i i i • i t*- 

• i i i i i i < > * 

i i i r . t i i ■ i 


CSA n- WP n’ FU n- AHT n- NHL n' MFF n- FFC n- F1 n’ F2 n> 


FLL *WP n 
n n 

4.0*CSA n 


where 


CSA- 

WP i 

FLL i 

AHT i 


NHL i 


MFF 


i 


FCC i 


F1 i 


is the fluid flow cross sectional area for fluid lump type i 
is the fluid wetted perimeter for fluid lump type i 
is the fluid lump length for fluid lump type i 
is the area for heat transfer for connection for fluid lump 
type i (usually WP*FLL) 

is the number of head losses for fluid lump type i if input 
as a real constant 

is stored as the relative location in the array data for the 
user input array of head losses vs Reynolds number if input 
as AXX where XX is the array number 

is the code to determine the method used for calculating 
friction factor for type i. If MFF = 0, the internal methods 
are used to calculate friction factor. If MFF = 

AXX, XX is an array (the relative location is stored) 
of the Friction Factor vs Reynolds number, 
is a constant to be multiplied times the friction factor 
for type i 

is a code to determine the method for calculating convection 
heat transfer coefficient for type i. If FI is real, the 
internal equation for flow in a tube is used and FI is the 
laminar fully developed coefficient. F2 is the laminar entry 


length coefficient. 

If FI = 1, F2 is AXX, XX is an arrayfstored as the relative 
array location) of Stanton Number vs Reynolds number array. 
If FI = 2, F2 is AXX, XX is an array (stored as the relative 
array location) of an array giving heat transfer coefficient 
vs tube flowrate. 

F2. is described under FI above 

i 


The tenth item in the list for each type is FLL*WP/(4.0*CSA) which is the L/D 
for the type. This item must be calculated and stored for each type during 
the preprocessor phase. 

The TYPE array is shown in Table B-I for the sample problem. 


4.0 Other Arrays 

Eight arrays must be set up for the flow problem in addition to the 
three primary flow problem arrays discussed in Sections 1,0, 2.0, 3.0. 

These arrays are each in a separate labled common block to provide ready 
access to them for user input and output in the user logic block. The 
labeled common block and the array name for each is given below: 


/WD0T/W(LT) 

/PRESS/P(LP) 

/FL0WG/GF(LT) 

/VALVP/ VP ( LV ) 

/WDOTI/WI (LP) 

/FL0WR/AFR(LT) 

/DELTAP/DP(LT) 


Array of flowrates per tube 

Array of pressures per pressure node 

Array of flow conductors per tube 

Array of valve positions per valve 

Array of imposed flowrate per pressure node 

Array of added flow resistance per tube 

Array of pressure drops per tube 


/FDIMNS/NTYPE,NSYS,NTB,NP,NV,NFD 

flow problem 


- Dimensions for the 


The dimensions in the above arrays are as follows: 
LT is the largest tube number 
LP is the largest pressure node number 
LV is the largest valve number 
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TABLE ft- 1 FLUID TYPE ARRAY 


CSA 

WP 

FLL 

AHT 

NHL 

0.001008 

0.1125 

12.0 

1.35 

WIM- 

0.000938 

0.36 

3.25 

1.17 

117.0 

0.001008 

0.1125 

5.0 

.5625 

0.0 

0.853E-4 

0.0328 

0.25 

.0082 

2.49 

0.001008 

0.1125 

20.0 

2.25 

0.0 

0.001008 

0.1125 

2.5 

.281 


0.001008 

0.1125 

50.0 

5.62 

0.0 

0.001008 

0.1125 

7.0 

.7875 


0.001008 

0.1125 

2.0 

.225 


0.853E-4 

0.328 

0.25 

0.0 



FFC 

FI 

F2 

FLL/D* 


1.0 

1.0 

1.0 

334.821 

10 

1.0 

1.0 

1.0 

311.834 

20 

1.0 

1.0 

1.0 

139.509 

30 

1.0 

1.0 

1.0 

.24E-6 

40 

1.0 

1.0 

1.0 

558.036 

50 

1.0 

1.0 

1.0 

69.75 

60 

1.0 

1.0 

1.0 

1395.09 

70 

1.0 

1.0 

1.0 

195.31 

80 

1.0 

1.0 

1.0 

55.80 

90 

1.0 

1.0 

1.0 

L__ 

.24E-6 

100 


D = 4.0 X CSA/WP 















The variables in the FDIMNS labeled common array indicate the size of 
various aspects of the total flow problem. The following values are assigned: 

NTYPE - Number of types 

NSYS - Number of systems 

NTB - Number of tubes 

NP - Number of pressure nodes 

NV - Number of valves 

NFD - Number of spaces in the FLOW array 
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APPENDIX C 

USERS DESCRIPTION FOR PLOT PROGRAM 


This Appendix presents user descriptions for a SINDA plotting routine, 
FLOPLT and a tape combining routine, MCOMB. Both routines are available on the 
ES3*SINDA Secure File. A brief description of the routines and the user input 
description is given below. 

FLOPLT DESCRIPTION 

The plot routine which is available on *SINDA can be used with a history 
file from a previous SINDA run to generate microfilm output. The items available 
for plotting are (1) pressure drop for each tube, (2) pressure for each pressure 
node, (3) valve positions for each valve, (4) flow rates for each tube, and (5) 
temperatures for each temperature lump. Each of these items may be plotted as a 
function of mission time. The user specifies the grid time range to be plotted, 
a time label, and the itmes to be plotted. A number of history files may be com- 
bined prior to plotting the results. The user has the option of averaging any 
portion of the plotted curve and of specifying the range of the ordinate axis. 

The system control cards and the data input card for FLOPLT are 
described below: 

SYSTEM CONTROL CARDS FOR FLOPLT 
0 RUN 

§ QUAL ES3 
@ ASG,A *SINDA 

0 USE 7, XXX (First file to combined) 

@ USE 8, XXX (Second file to be combined) 

Add additional USE cards as required for files to be combined. 

@ MAP *SINDA, S INFLOP LOT/MAP , TPF$ . RUN 
0 XQT RUN 
Data cards 
0 FIN 
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FLOPLT DATA CARDS 


Columns Format Title 

Card 1 {Title Card) 

1-72 12A6 TITLEA 


Card 2 (Parameter Card) 


1-10 

F10.0 

TA 

11-20 

F10.0 

TZ 

21-30 

FI 0.0 

TPG 

31-35 

15 

ITMX 


36-40 15 MPNT 


41-45 15 NTP 


46-50 15 KT 


51-55 15 INC 


Description 


Any 72 alphammeric characters to be used 
as heading for each frame of plots 

First value of time to be plotted (hours). 
Last value of time to be plotted (hours). 

Time range for each grid. Number of grids 
drawn will be (TZ-TA)/TPG. (If TPG is 
left blank, the job will terminate.) 

Time scale Table: 

= 1, "SECONDS" 

= 2, "MINUTES" 

= 3, "HOURS" 

Any other value, "*****'' 

Print control code 

= 1, prints information to be plotted 
while loading the plot tape 
f 1, will not print information to be plotted 

Number of tapes to be combined. Use a 
negative number if start and/or stop times 
are specified on Card 3 for any tape to 
be combined. / 

File number to which file to be 
plotted is assigned. If left blank, file 
23 is assumed. The combined file is 
assigned to this unit. 

= 1, every time point and associated data 
value from the tanes tn he combined 
will be transferred to the combined tape. 

= 2, every second time point and associated 
data values will be transferred to the 
combined tape. 

etc. 
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Col umns 

Format 

Title 

Description 

56-60 

15 

IUNIT 

Logical unit number to which first tape 
to be combined is assigned. If left 
blank, unit 7 is assumed. 

61-70 

FI 0.0 

ASTRT 

Beginning time for averages (hours). 

71-80 

F10.0 

AST0P 

Ending time for averages (hours). 

Card 3 

(Required only if NTP < 0. 

See Card 2 columns 41-45) 

1-5 

F5.3 

XSTART 

First time point from first tape to be 
combined which will be transferred to 
the combined tape. 

6-10 

F5.3 

XST0P 

Last time point from first tape to be 


combined which will be transferred to 
the combined tape. 

Repeat XSTART and XST0P in five column fields for each tape to be combined. 
Card 4 (Item Card) 

1-5 15 ITEM The item number to be plotted. Use a 

negative value if this item is to start a 

Repeat Card 4 for each item to be plotted. 

Card 5 

1-80 Blank 

Card 6 

1-80 Blank 

If additional history tapes are to be plotted, repeat Card 1 and subsequent 
cards for each additional history tape. 
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COMBINE ROUTINE DESCRIPTION 

The combine routine* MCOMB, can be used to combine history files 
into one history file prior to its being plotted or being compared to another 
file. The combined file which is generated can be saved for future use if required. 

The user selects the frequency with which the time points and associated data values 
on the original files are added to the new file. That is, every time point on the 
original file can be added to the new file or every second, third, etc., point can 
be added depending on the requirements for the combined file. 

The conioine routine is a very useful feature if several history files 
are generated on a long mission run. By combining these files before plotting, 
a continuous plot of the mission can be obtained. The convenience of the combine 
routine can also be observed when mission runs made with different time increments 
are compared. Obviously, the run made with the smaller time increment will take 
more computer time than the run made with the larger time increment, and will probably 
require at least one 11 res tart". In such a situation, there would be two history 
files with the smaller time increment to compare to one with the larger time increment. 
The two files with the smaller time increment can be combined and then compared to 
the file with the larger time increment on the same run. 

The system control cards and the data input cards for MCOMB are described 

below: 

SYSTEM CONTROL CARDS FOR MCOMB ROUTINE 
@ RUN 

@ ASG, A ES3*SINDA 

0 USE 7, XXX (First file to be combined} 

@ USE 8, XXX (Second file to be combined) 

Add additional USE cards as required for files to be combined. 

@ MAP ES3*SINDA, MCQMB/MAP,TPF$.RUN 
@ XQT RUN 
Data cards 
@ FIN 


C-4 


Columns 


Format 


Title 


Description 


Repeat XSTART and XSTOP in five columns fields for each file to be combined. 
Card 3 (Required only if K0DE2 > 0. See Card 1 columns 1 6-20 

1-10 FI 0.0 ADD Time to be added to each time read from 

first file to be combined. 

Repeat ADD in 10 column fields for each file to be combined. 


